“Calhoun 


Institutional Archive of the Naval Postgraduate School 





Calhoun: The NPS Institutional Archive 
DSpace Repository 


Theses and Dissertations 1. Thesis and Dissertation Collection, all items 


1998 


A Pareto frontier for full stern submarines via 
genetic algorithm. 


Thomas, Mark W. 


Monterey California. Naval Postgraduate School 
http://hdl.handle.net/10945/8803 


This publication is a work of the U.S. Government as defined in Title 17, United 
States Code, Section 101. Copyright protection is not available for this work in the 
United States. 


Downloaded from NPS Archive: Calhoun 


Calhoun is the Naval Postgraduate School's public access digital repository for 
| (8 D U DLEY research materiak and institutional publications created by the NPS community. 
a Ser Calhoun is named for Professor of Mathematics Guy K. Calhoun, NPS's first 


NY KNOX appointed — and published -- scholarly author. 

>. LIBRARY Dudley Knox Library / Naval Postgraduate School 

411 Dyer Road / 1 University Circle 
Monterey, California USA 93943 





http://www.nps.edu/library 















O | + 
a MDE Isia 
tw 


wet 












ae 
ite) rom OL 


Us 
















pi TRT - 











































































































































































































































































































































































































ERA A A 
% te a KAES Para 
A L 
CEET | A MA be i 
Dia CATE Bet KN? nn SR is 
ro Suton per ge KA He Pr CORRA by Pe Bee REA SECETEI ) E E 
AS] PAPUE E EE TOMS Pi Pat ye aa for A re A PNE 
E E DATA ARA a KRR t. ee bay i EAN va k 
2 PE] Fe of i Lp 1M ot ae PE i Ano d Ptr i | p 
‘ 7 k ALERTA O ES Ue th) 4 cy G] REP AO A 
t FOE Hor 4 Pr a q Pred t 
d PLATOS A Ya YM 0 hoy z ay re eee ty 
A ACA A eri »h 
; atas a A ATAN waa IER A ON yr Pan 
ig PY loan, Mo. ne Pe re my he ET ARO 
o RANA EPA Med NE Pe Pn a 
a Stare rae ar Oh aa As by 
E , a b i EL 7 Corer aay ui 
A! AA EPT 
TOFT det IAN praia ad at) 
; E PA af STO Ol 
ar RO ea ARRI Y dm 
KUTETA PRN AA 
; A E aos AT 
l m A l ORTA A fea slan E ETA: 
$ i Re e ati Ce TIN 4 Pik rae fal 
, | APEPI N 
A Ñ anto LI UPA E CA Li a 
“ y A] TNH 1 TARRI 
A + eye Vin, 
0 A 
TEn .a 
a Y de 
: y A 
A 7 A 1 ed, ON mu Le 
; A S Merits One eee era UN j 
Chat} o 9 0] f A Spee a AMADA 
A TIA ye E 
a A E , y ee oP sa, AAA A Po AT oa Es Ed 
Pa TA, TREENERIT Soia 8s! fh SUPA ety ty ket hy 
A LON e e an KTERE De a E Panir » 
an a , n a AA E Hi A SUR? iad Fe AT oe 
He AI RA e n S GART ILES] HE LEER TEN hst PAY a 
A ; k; AY AA ATA SR A 
oh eres J PEREA NO TA TA De A 
i e P Ah o’ - ADT AT IAN 
Ar, AIN AAA 2 
ee Pan p AN 4 A IN A CIA TE AY] hor phe“ LS a ARTA 253 
A : CE A AA EA A CRI ie AAAA 
ne j AN AOS Por re Cores ee AN eha, BnS 
4 g A PL ow EE TOYE ACTELE ES DN ` 
; nie Au , b : 3 i i, 
n o h ALLT 
; CL ICR RTRT IERP E PEETA Sorte h 
rek y na Chics ati oe oe dad vegan PAE A TERA EITI 
yal EM A EFEKT IPTA ES MEA H 
E ; aie daa ay 1 A RET a OD A A 
A A A A ee CEEE narar G AAN 5 TA 
De DA ea t TE ETE ASAS 
A y da ea bos , ery i 
A al EI 1%) Maia Vallas 
ta A 
0 
t s 
© à 
LCi 5 
.. 4 ef 
aio e 





LS 
Ash a EARN ER 





















































































































































































































































































































y ts OO 
A E acd ion TA 
TES E ETY gd sd 
Van a EEE EN PE Peres ' 
te Minga y Milpa. AE AS 
T ATI ey | AAN 
A Ya ` | T een 894 ie, neta TY A AA 
ec ï A OSA DER KN A MD 
s a ry $ Vds ino AN ey 
A on 1 en vr AI O KPE O ” 
, A > Ar Vel cn UA PLTA FR 
A POE! LON z h A T S be 3 Paes) EIA et 3 
9 A coe wth ote ga Naty dt Is | ES Ue O STE) EA Pid Pave Ot, PLP De ayo ee 
o sus ste eee TE CN IMA PAR NS IS A eye, 
y sos o AE $ O tng es IA A IA Mba IR e 
bye Sime lady ‘a AI A e lsrs e, AA Ng n X KORCARI] edo ca ia de 
s E , Be HA OH JR O + O O AI ee Peewee PE ATA HITOS AIR 
RNE SS T O CA ie ERTZ E TR MA Ma Sul UTERET SPT PE BrE 
ARA O AA E ae Aid CW Fire ree Tara PC Pe epee ANS ds 
an f > ST TT A PA AN Ly 
$ ae) IF NHS nla AT bet teen er eee 
pe as NEL RAE TA a Ae Jie 
k AN tapo? SIE EPA TI La o or 
aa ds tet, fel ih Ha ka 
ES i os Cr PAE ray 
ITA TÍA pi ‘ A! 
"eh ae b: PESA + K Bda rh 2 f EES P 
E A e CRUE ODE E Ss e 
an y eS ow ah A re Ok ee X Wadia a UEM SR 4 A TEA 
3 A A A y TAN EA PR a rad Io 3 A DT Y 
A Pe) » CA PA WV A O tN TEN Me A Pak tig eras ee ee ply 
d SU N A: Mah stone 1Y ETA A AT PET tye, TASTE MA y 
i Ag 1 Mis 4 ki- A y ier h tgi E EST PEIROT TYE a pa > rias 
A IA wos e { Tega AA PS ESE se Peep te peel a LE 
a y t) w ELR EI A DETT dé ts PTA LE 
3 TETIDO wr OIA eA A EE 
e TT AS . “i PFGE oe a A A EI A 
= Tete seme ye i Lora A TAI cri ee, ath: Li 7 E RS DLA % 
‘ A o E ae Jac Cd A G G MTY ees TT A a J Fe an ec a 
« i A o II T s das ja > A A Ki SLETT EES , n >= 24 H Asan p E teow YW a CERI 
y y IS a IS O fey i tat + ' CS ` 1 (3 4 TES À 7 A ER TO Xx Pads A APEA ds * par A 
A A CIRO OE TO CI at eT | Lhe ee Loan q. ASS » E pd pa 
A oo. y O be W» Las » 4 G b Lipa joa | wd e .. y e t te 7 k tst iseda ng $ 7 
‘ aes | Sy. » 4 .. . ¡Ba MAD | .na y O ot Do ' A Ame srl Mn CAE TER 
* a a .a (PR A A E mn n % O > TS (E "y t bie at A j PEA 
0 r Dun g y A AT 3 0 ? RY Por] A eb. Aa D 
n 7 ae Fy feai va La ae Oa. 4 
ÓN i ATA AS 
AP y Ona IA A PS hs Pure y 
CA IIA II A nat Lea M4 S 
e A O A A mS € 4 
1è iye’ + tw 194, i * ta 
Ys tree 
A PPT Soe Newey hy 
rar A AA sebit o 1 A à 
ETET a DI ʻi E 
CA re e. 4 
T ct S OA MAI : 
o mta PAE AREL s 
A y TA A z ETE e3 ee 
TA An EAN 
AS IT O r 
E p ba ” hh 
A A y Ai vo 
te >. 19 Var ona m.s 
0 CEE A T "is abr o. 
$ °- ere pate . ay 
A wi CERT een à H 
À E OS CI ERA 3 
i. s T ES E ICA Te po TS 
f tr mu d EE IE TE AS t 
ML) Yat ON A T A 
A A p 
` 
D 
P 
A 
D 
o 
g 
AG: 
7 
An 
. 
‘ 
A 
` 
o 
A 
y 








hes fare 
Petts ie * 
CS Ye a 








mm. 2 









Ly 
v 
13 y 
kl Be] 
R Atte 
A w'¢ 
ae 
¥ 












« Ca Dei 
OE 
e AT 

PURO 






bea IEA 
cil ee 
P HAARA suig 
E i “f° EET IA E @ 
Peso a ron, 
ISS s FR CEE BY ea 
ae tek A re E bg al ir “TM pag it 
y O A AS “ 









































































Mote Mg be 
ee IN y PA o ae! aupa E y 
LEAL Es AE ar E a PEE PELIT PRIFT ee ae 
(OPET) i a gtp A A IN k nN prats bra y a wg F» PT EE AAEE juig 
CEE ES SC EET LEN J Gig ote gw e , pa 

n A H 


A A To 





A O Or ae po 
>] io y 
t 


ROA TS 
[sy aa ‘ =e 
EE 


E 
z E A 
AR EA E 2 AA CS AT 
ET AT e A E eld 
HO ha AS rum ME gee ae UB a Ie cad 
dad, AD eg RO 
E AT AR LIFE RI Sf JA AS ETE bee Gnd 0 ety ut ag hy shag e RTE 
4 ITA rs) a erro) SA y e AA y 
A wan Oa Gee T ITET s j LELT N E] AI ISS EN 
NoE rp th ioe bid epa ‘pty 









CIT | tee tr f 
awe Aane an A Ur, 


























de! Cots 

[RIA ES Y A 

IA A A A PRATS FS ogres vee 

ee ae AAA Wee ate a igs > PRTA 

y A rat SOL HE A LELE CSEL ENT LE A ATAN 

A? tq SA ad ce Pay ert “vas gtd Nee Depa 
J Ae z taps terjun ITA ee HORUT 

yo. ee o DEn E n yk? arid EAN ud 
RT A AS o a Zeas SES t AE 
Poe tuya $ ELE 










oe) PoP oo da 
cl AAN TAS qe 
re Leas Be) y aka oe (ead AE dis 
rr EMAN Pp wy eee 
AS Ra de 
MAA KERREN SPs ieee nar 
ar vajej EEU ná 
E TA AS 


> 
LS or ihe erp 
Ce 


qu E 














E 
Ai br a 
ee he ek E 
RITA ad PE Ae le 








5 eae “a UA eies y a i 



























DALS 
e CNA AA CE Ds a “ete at dete PICS 
g GA uf y ee € 
A E S ae K4 ai 

ig 4 





E Eb bs Ms it EA ae Pu 
1 Fe a yh A E p ex Cait 
ed'o y A oo E 


e [o 
ie NER i A t if ahn ee 
». A Ls nadas yl ay 





en Sr Pee sé 
q 




















































































su TA Pot ado 
ORO hy tars RA Fa ed a a PERAS pá 

A ASS ay Abo eye, EA AF 
y CN 5 M E FHA VERE LI yn 

aa A! A Y UT ION a da 
4 ue STA ld id at hd bat J =p 
r T ee Me y NU 
E OA | Mas 1 he» ER CAEN qe 
7 a J TA E ey rae fa bee es ak faari, mM 
"y e E par E g tn Sy È ITET S TE hy Pa ng O x rie CAS Aapan 
.t 1 e? f Nie REO AA IO DIE A AA T ` EDT IPD NE Pe k 
EES n EET a LITA ryo ATE IA A St ET 
er ee aE oJ Pa vey 1 AAA „>! 
e yh pa 


py Pia ne NE G 
RRT EIP 










Orr Pca 


A O 
2 Loi | Ds afl 


rá ARTIC 
> tt sae gt 
















ir, ae A y; KIS red 
MUI A RIEA PA TE 
US A 
HRU t 
efor 









Ca SAE, Tee or, 
| e y rey 
car Ls 


E A S 
naa rr HA 

ee ety a 3 ra 
MAT RERPI ahi 


ot 



























Ree E 3, OR y4 
M'i LEDT A Jw rp. za 
r SrA CKE E rhb a] rh) a Le ha 
mata zs TN MT PA Ca 
07 ar CTA] POD AS EA a ot 
Nth a nn hiy yhaa PALO Ray 
DATEI vara ory Da ATA ES 
E ro meh dieta! 10d, iy AAA MT 
ALH ELY T eat 
e 


A CH et et are 
pan gav 
















Y 
E? 

















qe 

Preto a pene 

VENA le Ay yer E y AY PE 
CIA] 


ES 




























mus 
ADA aha Li 

































TES AS 
g TORET TE 
7 E NEO ey REAT Fi pt, ATL oe 
Ag PG 8, yank Pa Sar Ds Y, E SES A CAT 
a Le eps N e Pit ree y gr A a ARAS e 
Ef Abe aes ong DOTA thy 
i a § Crh eas 


a 
oO ad 
TE 








3a 
LEE PEITTI 
PN e a 
een Matt ee tap 


ap ay 


ed ae KIHE Ry 
SAS pa 

Se 1 Pat. oot 
E è T $ 01 
va a pp A 



















AAA 
RR. A a, RES a.y TORE at 
ty PU re CA wats: vate Le ay ee 
E ER Ape Teed! aco nag ett er 4 atta PA ter bas NR fs 
f va TI T UP 
tes try mi nu vir 


ES rE 
AE 


ZB pds asd so ott a 
PE 


DUDLEY KNOX LIBRARY 
NAVAL POSTGRADUATE SCHOO! 
MONTEREY CA 93943-5104 


MUDLEY KNOX LIBRARY 
A. POSTGRADUATE SCHOOL 
SONTEREY CA 93943-5101 








A Pareto Frontier for Full Stern Submarines 


via Genetic Algorithm 
by 
Mark W. Thomas 


B.S. Electrical Engineering, Oklahoma State University, 1984 
S.M. Electrical Engineering and Computer Science, 
Massachusetts Institute of Technology, 1996 
Naval Engineer, Massachusetts Institute of Technology, 1996 


Submitted to the Department of Ocean Engineering 
in partial fulfillment of the requirements for the degree of 


Doctor of Philosophy in Hydrodynamics 
at the 
MASSACHUSETTS INSTITUTE OF TECHNOLOGY 
June 1998 
(¿Mark W. Thomas, 1998. All rights reserved. 


A FP AVON VS 
UB — PS 
à 128 


the WAC S YV , æ 





A Pareto Frontier for Full Stern Sabin KNOX LIBRARY 
via Genetic Algorithtif vir OSTGRADUATE SCHOOL 


he EREY CA 93943-5101 
Mark W. Thomas 


Submitted to the Department of Ocean Engineering 
on May 12, 1998, in partial fulfillment of the 
requirements for the degree of 
Doctor of Philosophy in Hydrodynamics 


Abstract 


An exploratory design and analysis code for underwater vehicles with ducted, multi-stage propulsion 
systems is developed from an existing version. Boundary layer modeling is added and used to 
predict flow separation as well as to estimate the effect of wake fraction on propulsive efficiency. 
Force balance (i.e., convergence to a self-propelled condition) is achieved by automatic variation of 
advance coefficient. The force and boundary layer calculations of the revised code are validated 
using published experimental data. 

A slightly modified version of the code is then used as the evaluator for an original Pareto genetic 
algorithm. The algorithm seeks the code’s Pareto (non-dominated) frontier in terms of usable hull 
volume and propulsive efficiency, with the intention of investigating the feasibility of so-called “full 
stern” submarines. Three different methods of Pareto selection, drawn from the current literature, 
are installed in the algorithm and compared in terms of their ability to locate and define the frontier. 
A new concept in evolutionary computation—non-interbreeding competitive species—is introduced, 
allowing simultaneous optimization of four incompatible propulsor configurations. This produces a 
feasibility frontier with optimal propulsor configuration as a function of stern fullness. The ability 
of the algorithm to define a three-objective Pareto surface is demonstrated, by including minimal 
cavitation as a third objective. 

The results provide evidence for the viability of full stern submarines, demonstrate the utility 
of genetic algorithms in obtaining Pareto design frontiers, and show that Pareto optimization is 
preferred to scalarized multi-objective optimization in general decision-making. 
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H* -- boundary layer kinetic energy shape parameter 
y e traditional advance coefficient 
JB = alternative advance coefficient 
ks empirical constant for cavitation correlation 
<T Sr traditional thrust coefficient 
“go ope traditional torque coefficient 


length of a bit string in number of bits 


vehicle length 


15 


m number of times a given schema appears in a population 
n propeller rotation rate (==) 

O order of a schema (number of non-“#” characters) 

p pressure 

Po probability of crossover 

Dm probability of mutation 

es total far-field fluid pressure 

P population size 

P propeller pitch! 

P/D propeller pitch-to-diameter ratio 





shaft torque 


i radial coordinate 

Rpg maximum vehicle radius 

Re LL Reynolds number 

Re. Ha blade tip chord relative Reynolds number 

T thrust, in general 

U tip speed (used in cavitation correlation) 

U velocity vector 

Ue magnitude of boundary layer edge velocity 

Uj iP component of velocity 

Ur circumferential mean self-induced tangential velocity 
u local self-induced tangential velocity 

Ua total local induced axial velocity 

Ut total local induced tangential velocity 

V, reference velocity (vehicle speed) 

W; relative fluid velocity at blade tip 

Wo edge velocity of end-wall boundary layer 

T longitudinal coordinate (Huang experiments only) 

í boundary layer curvilinear coordinate parallel to local surface 
Y boundary layer curvilinear coordinate normal to local surface 
ž longitudinal coordinate 

Z number of propeller blades 

e circulation per unit length (4) 

iE irt circulation on a lifting line ($) 


l Propeller pitch is the distance the propeller would travel through its medium during the course of one revolution 
with no slip. It is generally a function of radius, and may also be defined as 2rrr tan $, where ¢, the pitch angle, is 
the angle between the intersection of the chord line of the section and a plane normal to the propeller axis. 
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pu dr. general circulation (E) 


length of a schema in number of bits 


7 iene — uz) dy boundary layer displacement thickness 


DSB D A 


* 


SD 


IS nN a a O; 


ee 


value of y at which lu¿(y)| = 0.99u. 
tip gap to blade chord length ratio 
tangential coordinate 


i} a —uz)uzdy boundary layer momentum thickness 


0 

Jy (u2 — uz)uzdg boundary layer kinetic energy thickness 
vector of objective function values 
Goldstein reduction factor 
air content correction factor 
ratio of tip gap to maximum blade thickness at tip 
kinematic viscosity ES 
fluid density 

Pmin cavitation index 


viscous shear 
viscous shear at hull surface (wall) 


vector of design parameters 
=) 


hull volume, less gear and shafting estimates (see Appendix B) 


angular velocity ( 


wild card character in a schema 
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Chapter 1 


Introduction 


This research is undertaken to support two independent theses: 


1. The volume and arrangement advantages offered by full stern submarines may be 


realized, without efficiency penalties, by properly designed ducted propulsors. 


2. Valid design optimization in the presence of multiple objectives must be based on 
knowledge of the realizable solution space. This requires definition of the Pareto 


frontier, which may be obtained by genetic algorithm. 


An existing exploratory design and analysis code for submarines with ducted propul- 
sors is revised and upgraded to perform the calculations needed for optimization. 
This code is then used as the evaluator in an original Pareto genetic algorithm, with 
dual objectives of propulsive efficiency and usable hull volume. The algorithm seeks 
the Pareto frontier by constantly pressuring a random initial population of variants 


toward improvement in terms of both objectives simultaneously. 


1.1 Ducted Propeller Systems 


The potential benefits of enclosing a marine propeller in a duct were proposed as 
early as the 1930’s [58]. Depending on the duct’s shape and construction, these 
benefits can include greater efficiency, mechanical protection or reduced noise. In 
some cases, however, the additional design, production and maintenance costs of early 
applications offset any benefit obtained. This concern remains relevant to some extent 
today, and is why the majority of marine vessels continue to use open propellers. 
One rather obvious benefit of a duct is that it can prevent propeller damage and 
fouling when the vehicle is operated near obstacles, in shallow water or under ice. 
When protection is the primary function, a duct is perhaps more accurately called a 


shroud, although terminology varies in the literature. Protective shrouds are usually 
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Figure 1-1: A ducted propeller on a surface vessel. Taken from Gillmer and Johnson, Introduction 
to Naval Architecture, Naval Institute Press, 1982. 


thin and straight in cross-section and are not intended to have significant hydrody- 
namic properties. They are aligned more or less with the direction of flow—so as to 
generate no force—and may be slotted or perforated to allow through-flow. These 
are fairly common on small submersibles such as autonomous underwater vehicles 
(AUVs) and torpedos. 

A duct may also be used to decrease the mini- 
mum propeller diameter for a given thrust; that is, 
to increase the maximum load that the propeller 
blades can carry. This was the primary motiva- 
tion for early applications on tugboats. An open 
propeller is normally designed to be lightly loaded 
at the blade tips; any significant pressure difference 


across the blade at the tip will cause the flow to 





spill over, resulting in a tip vortex as shown in Fig- 
Figure 1-2: Tip vortex cavitation on 


an open propeller ure 1-2. The production of a tip vortex requires 


energy and therefore causes a loss of thrust as well 
as possible cavitation, noise, vibration and erosion [10]. However, if a duct surrounds 
the propeller and the clearance between the blade tips and the duct inner surface is 
small, a pressure differential can be maintained across the blade tip. The result is 


an increase in the load-carrying capacity of the blade, a reduction in the required 
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propeller diameter, or some combination of the two. 

Both propeller protection and increased load capacity are possible with simple 
shroud-like ducts. Further benefits are possible if the duct is given a thickness and 
camber distribution and/or is operated at a non-zero angle of attack, in which case it 


becomes an annular hydrofoil (see Figure 1-3). The duct is then capable of altering 


chord (c) 


thickness distribution 
mean camber line 


angle of attack 
(a) 


camber distribution 


F(x) 





Figure 1-3: Typical duct cross-section with parameter definitions. 


the propeller inflow velocity as well as generating lift.1 This is illustrated in Figure 1- 
4a, which shows the stern of an axisymmetric vehicle operating at constant speed. 
Thrust provided by the rotor exactly balances the viscous and pressure drags acting 
on the vehicle. Streamlines generally follow the contour of the hull, perhaps with 
some small contraction imposed by the suction of the rotor.’ Figure 1-4b shows the 
same configuration with the propeller now surrounded by a thin duct, represented by 
its two-dimensional cross-section. If the duct’s mean camber line lies along a former 
streamline as shown, the duct generates no lift; its only effect on propulsive efficiency 
is some small axial force due to its viscous drag. Now let the duct rotate slightly 
about the propeller tip such that it attains a small positive angle of attack relative 
to the inflow, as in Figure 1-4c. The Kutta condition dictates negative (clockwise) 
rotation in the surrounding velocity field to prevent infinite velocity around the sharp 
trailing edge. This rotation may be quantified by a contour integration of velocity 
around the duct; it has dimension of length squared per unit time and is known as 


lIn very general terms, ducted propellers which decelerate the inflow are quieter and less efficient than equivalent 
open systems and are referred to as pumpjets. Ducted systems which accelerate the inflow are noisier but more 
efficient and are referred to as nozzles, or sometimes as Kort nozzles after their inventor, 

2Streamlines aft of the rotor actually follow a helical path; those in the figure may be thought of as a centerline 
slice through the wake. 
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rotor blade 





Figure 1-4: Effect of duct on propulsion system. (a) Non-ducted system. (b) Thin duct aligned with 
former streamline. (c) Duct loaded by placing it at a non-zero angle of attack. 


the circulation ([) of the duct. Application of the Kutta-Joukowski law 
Lift = pU xT (1.1) 


reveals that the duct is now subject to a lift force normal to the direction of the local 
flow. Viscous drag is still present and acts parallel to the camber line, along with 
additional induced drag due to the lift [54]. The radial component of the resultant 
force has no effect on the vehicle’s propulsion system, since the duct is axisymmetric, 
but the axial component F, must be compensated by rotor thrust in order for forces 
on the vehicle to balance (note that if the angle of attack were negative, the lift would 


be in the opposite direction shown and the duct would provide thrust). 
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The effect of the loaded duct on propeller inflow is indicated by the altered path 
of the middle streamline; this alteration will affect propeller forces and thus overall 
efficiency of the vehicle. The blades, like the duct, are lifting surfaces and are subject 
to the Kutta-Joukowski force. T'he difference is that the vorticity vector of the rotor 
blade shown in the figure lies in a z-r plane, while that of the duct lies in a 2-9 
plane (the blade I, if shown, would be edge-on in the figure). According to Equa- 
tion (1.1), any variation in the axial (z) component of velocity at the rotor alters the 
torque required to turn the shaft, while variation in tangential (@) velocity affects the 
thrust produced.’ Both of these quantities affect propulsive efficiency. A loaded duct, 
therefore, affects a vehicle’s overall propulsion characteristics in two primary ways: 
indirectly, through alteration of the propeller inflow, and directly, by generating lift.. 

Whether and how the combination of these effects can be made beneficial to the 
vehicle is not obvious, because the mechanisms described above are coupled with each 
other and with body forces as well. The increased velocity due to the suction effect 
of the propeller will lower the fluid pressure at the stern and increase the pressure 
drag on the body; this added drag may be augmented or diminished by the duct 
circulation. The propeller and duct loads, however, are dependent on the inflow 
velocities presented to them, and inflow direction is very much a function of the body 
shape at the stern. Stern fullness—the abruptness of transition from maximum body 
radius to the rotor hub radius—thus affects the velocities induced by the propeller 
and duct which in turn affect drag on the hull. It may be possible to manipulate this 
interaction among the hull, duct and propeller by adjusting design parameters. If 
so, it would be of interest to know which combination of parameters—that is, which 
propulsor configuration—requires the least power to propel a given hull form at a 
given steady speed. Such information would be particularly useful in investigating 


the feasibility of full stern submarines, which are discussed below.* 


1.2 Full Stern Submarines 


The full stern submarine is a conceptual departure from traditional submarine design, 
involving a relatively abrupt transition from maximum hull radius to propeller hub 
radius at the stern. Figure 1-5 compares a notional full stern profile with that of an 
ideal, minimal drag profile and the common parallel mid-body/tapered stern profile. 
The ideal profile represents a trade-off between friction drag and pressure, or form, 
~ 3 Assuming that the propeller continues to rotate at the same speed, of course. 


1Capt. Harry Jackson provided the author's first exposure to the concept of full stern submarines at an MIT 
Professional Summer course in 1995 [50]. 
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parallel mid-body with full stern 





Figure 1-5: Notional full stern submarine profile, as compared to the ideal (in terms of drag) profile 
and the common parallel mid-body/tapered stern profile. The internal walls represent possible 
pressure hull boundaries. 


drag on the hull. Friction drag is directly related to surface area and fluid velocity 
at the surface; pressure drag is a more complex function of viscosity and momentum 
recovery. It is sufficient to note here that pressure drag is minimal on a long, slender 
body with a slowly varying cross-section. Such a body will experience relatively high 
friction drag, however, as it requires more surface area than a shorter body of the same 
volume. An optimum length-to-diameter ratio should exist, then, which minimizes 
the sum of these drags. This ratio has been estimated at 6:1, although the total drag 
curve in the region of the minimum is quite flat [10]. 

There are several reasons why the ideal profile and optimal length-to-diameter ratio 
are not common among operational submarines. Given some required volume, these 
characteristics it may require an unacceptably large diameter or allow insufficient 
length to contain internal machinery and systems. The non-constant radius of the 
ideal hull increases production costs and presents difficulties in both construction and 
maintenance (e.g., dry-docking). It also makes the shape of the inner (pressure) hull 
problematic. If the pressure hull conforms to the outer hull, as in the upper profile of 
Figure 1-5, internal arrangement and deck layout are difficult and some pressurized 
volume may be unusable. If pressure hull walls are made parallel regardless of the 
outer hull curvature, extra volume is introduced between the two hulls. This increases 
the total volume to be propelled, possibly resulting in an overall efficiency decrease. 


Such considerations have driven the majority of submarine designs toward a parallel 
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mid-body profile, such as the middle profile of Figure 1-5. The advantages realized in 
cost and arrangeability are assumed to offset the increase in drag, which in any event 
is slight due to the flatness of the total drag curve. 

Parallel mid-body designs generally retain the long, tapering stern shape of the 
minimal drag profile. Tailcone half-angles do not normally exceed 20 degrees [88]. 
Construction complications associated with non-constant radius are therefore not 
eliminated entirely by the parallel mid-body plan; they are simply avoided along 
most of the length. Also, the attachment of the long stern section to the pressure hull 
remains a structural challenge, increasing in difficulty as the tailcone angle decreases. 

A short, full stern section would reduce or eliminate these concerns. Also, since 
propulsion systems and other high-volume machinery are usually located aft, a full 
stern would allow greater flexibility in arrangement. The additional volume might 
also be used to increase accessibility for maintenance and/or to reduce the overall 
length of the vehicle. Reduced length would bring most parallel mid-body designs 
closer to the optimal length-to-diameter ratio. Such shortening has been investigated 
by Warren [88] and shown to be feasible, possibly even resulting in greater maxi- 
mum speed or propulsive coefficient.’ Full sterns also produce relatively large wake 
fractions, or “viscous shadows.” Propellers operating in these regions of reduced ve- 
locity may realize an efficiency advantage over an equivalent open-water system. This 
potential benefit is taken up in greater detail in Section 2.3.1. 

Despite the possible advantages of full sterns, they remain notably absent among 
current designs. One reason for this—the desire to limit form drag—has already been 
mentioned above. Another reason, more compelling, is to avoid flow separation. As 
the fluid near the body surface flows over the stern, it decelerates due to the increasing 
pressure.° If this deceleration reaches a critical value, the flow will detach from the hull 
and proceed more or less directly downstream; this condition is known as separation 
(see Figure 1-6). The likelihood of separation at a given speed depends primarily on 
body shape. Bluff, or blunt, sterns with rapidly changing cross-sections tend to cause 
separation where slender, tapered sterns do not. The consequences of separation 
are quite undesirable, and body profiles with the potential to cause separation are 
avoided. Separation causes a significant increase in drag due to the low pressure region 
it creates behind the vehicle. Fluid in a separated wake is relatively stagnant on a 
large scale, but is at reduced pressure due to the requirement of continuity with the 
Ma ande propulsive efiaicney, but is not relerred to asisuch becausejlt has 
no absolute upper limit. See Gillmer and Johnson [30] or Burcher and Rydill [10] for an overview. A more detailed 
analysis, relating efficiency to the various physical phenomena involved, is given by Dyne [23]. 


© Of course, it is actually the vehicle which is moving, but a stationary vehicle with fluid flowing past it is an 
equivalent situation and is consistent with most models in the literature. 
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moderate adverse pressure gradient 





Figure 1-6: Graphical comparison of attached and separated flows, showing how streamlines detach 
from the hull at the point of separation. Fluid in the separated region is at low pressure, causing an 
increase in drag. 


free-stream pressure outside the wake [67]. The propeller is forced to operate in this 
region of erratic local velocities, making efficiency unpredictable at best. Additionally, 
separation is generally accompanied by an increase in noise level. This, of course, is 
a particularly unwelcome effect for an operational submarine. 

Fortunately for proponents of fuller sterns, the stern shape is not the only factor 
aftecting the likelihood of separation. As previously mentioned, the suction effect 
of the propeller is felt upstream for a distance of a few propeller radii and a duct, 
if present, may also alter the flow field locally. It is conceivable, given the complex 
interaction among propulsor components, that ducted propulsors can be designed 
to delay or prevent separation which would otherwise occur on a full stern without 
degrading overall efficiency. 

This presents several interesting questions for submarine designers. There are 
benefits to be had by increasing the taper of the stern, but separation must be strictly 
avoided. Propeller and duct effects, which are interdependent and coupled with the 
shape of the stern, may delay or prevent separation. How abrupt can the stern 
transition be, and what is the best propulsor configuration for that profile? Propulsive 
efficiency is critical, particularly in non-nuclear submarines, where a twenty to thirty 
year life cycle means that operating costs will exceed those of design and construction. 


Does the maximum attainable efficiency decrease as the stern becomes fuller, and if 
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so, are the volume benefits worth the efficiency penalties? Ducts are known to affect 
propulsive efficiency and cavitation noise. Should a duct be used, and if so, should it 
accelerate or decelerate the propeller inflow? Guide vanes, or stators, can be placed 
forward of the rotor to provide efficiency-enhancing pre-swirl or aft of the rotor to 
recover wake energy in the form of tangential velocities. Any benefit so obtained, 
however, must be balanced against increased system complexity and friction drag. 
Should stators be used, and if so, should they be forward or aft of the rotor (or 
both)? 

The answers to these questions involve not only an unwieldy number of design 
variables but also two possibly conflicting objectives (volume vs. efficiency). This 
research, in addition to seeking answers, is also intended to demonstrate that such 
questions can and should be answered satisfactorily in any design process. Further, 
they can and should be answered without ad hoc or a priori assumptions regarding 
the relative importance of the objectives. These are general issues, transcending 
hydrodynamics and submarine design, and will be dealt with as such. A likely field 


in which to begin is that of multiple criteria analysis, discussed below. 


1.3 Multiple Criteria Analysis 


Decisions are not based on isolated, independent criteria. Any choice among alter- 
natives, be it the propulsor configuration for a submarine, foreign policy or what 
to have for breakfast involves prioritization of the consequences, either implicitly or 
explicitly. Certainly there are instances where a single criterion dominates (whether 
one should exit a burning building, say), but these are generally cases where no de- 
cision is required; the choice is obvious. In situations where a non-trivial decision ts 
required, its difficulty may be assumed directly related to the number of objectives 
(implications, consequences) involved and the degree to which they are conflicting. 
This assumption has a direct impact on engineering, or design, optimization. 

The impact is due to the fact that any optimization, including design, involves 
decision-making. In fact, a decision is usually required before the optimization pro- 
cess can begin; its purpose is to define the optimum. Take for example the linear and 
non-linear programming techniques which exist for the “optimization” (minimization, 
usually) of mathematical functions. Prior to invoking them, the solution must first 
be defined as the parameters which return the smallest value when evaluated by the 
function.’ Such a definition—though implicit and trivial in the case of mathematical 
function minimization—is more obscure in engineering optimization, where multiple 


“In a constrained problem, of course, the solution must also satisfy constraints. 
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conflicting objectives are likely to be present. The defining of optima in these situa- 
tions is complex enough to warrant its own field, that of multiple criteria analysis. 

Multiple criteria theory had its origins in the 1950’s; the earliest use of multiple 
objectives in mathematical programming is usually attributed to Kuhn and Tucker 
in a 1951 paper [57]. Interest in the field intensified in the 1970’s, as evidenced by 
the publishing of several new techniques and practical applications. As the number 
of techniques grew, basic procedural differences became apparent among them which 
allowed classification of the field. MacCrimmon [61], Hwang and Masud [49] and 
later Dlesk and Liebman [20] attempted to establish a general classification system 
by distinguishing between the terms multi-objective (describing an engineering or 
design-type problem, generally having at least a partially continuous solution space) 
and multi-attribute (describing a decision-type problem with a finite number of pre- 
determined alternatives), classifying both as sub-categories of multi-criteria analysis. 
Hwang and Masud [49] further classified then-current optimization techniques de- 
pending on the point in the process at which subjective preference information (the 
decision mentioned above, essentially) is required. This remains an important dis- 
tinction among techniques, and it is emphasized in a recent comprehensive review of 
the field by Miettinen [65]. The present research takes the position that it is advanta- 
geous to delay subjectivity as long as possible (it can in fact be avoided completely in 
rare cases). Prematurely applied subjectivity excludes regions of the solution space 
from consideration before the form of the solution space has been investigated. This 
can result in poor decisions, as will be seen shortly. 

Among the techniques falling under the general heading of multiple criteria analysis 
(listed here without regard to classification by subjectivity requirements) are global 
criterion, €-constraint, bounded objective, goal criterion, utility analysis, analytical 
hierarchy process and goal programming.® Each of these is a unique philosophy of 
what constitutes optima in the presence of multiple objectives. In choosing among 
these techniques (again, making a decision), one accepts a definition of optima. A 
definition combined with a search mechanism is an algorithm, which may be used to 
locate the defined optima given a mapping between parameter and objective space 
(1.e., given an objective function). 

Such algorithms require not only a definition of the solution but also a way of 
distinguishing between superior and inferior candidates during the search. Again, 
this distinction is trivial for single objective minimization—one solution is obviously 
better than another if its objective value is lower, and the optimum is the feasible so- 


8 Generally, whether a method is intended mainly for decision-making or for design optimization is not obvious 
from its name; both types are represented among these examples. 
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lution having the lowest objective value of all feasible solutions.? The global optimum 
is normally sought, although some problems require identification of both global and 
local optima (see Figure 1-7). 
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Figure 1-7: Comparison of optima for various types of objective functions. (a) Single objective with 
global and local minima. (b) Single objective subject to constraint. (c) Single discontinuous and 
disjoint objective. (d) Ambiguous situation with two objectives. 


Unfortunately, such simple definitions of “better” and “best” are not suitable in the 
presence of multiple objectives. The reason is shown graphically in Figure 1-7d. Two 
objective functions are present, Ff, and Fy. Although both of the candidate solutions 
(Y and S minimize one of the objectives, there is no rationale for declaring either 
solution “better” than the other.*” The ambiguity cannot be resolved by assuming 
the objectives to be independent and performing separate optimizations. Solutions 

° The definition also applies to single objective problems where the objective is to be maximized, by simply changing 
the sign of the objective function. 

10 Simply defining the better solution as the one which minimizes the sum of the objectives is equivalent to assigning 

them equal importance, which of course is subjective. Even if this were considered acceptable, there remains the 


problem of normalizing non-commensurate objective dimensions (the cooling requirements and weight of an electric 
motor, for example [84]). 


29 


Q and R, for example, both minimize objective Fy but neither is very attractive in 
terms of F). 

In fact, even the term “optimum” itself is ambiguous in the presence of multiple 
criteria, as evidenced by the number of published methods for defining it. This ambi- 
guity may be better appreciated with a simple example. Consider a team of engineers 
tasked with designing an optimal nuclear power plant, in terms of construction cost 
and safety, for a given power rating. They are hopefully not interested in the cheapest 
plant regardless of safety nor are they likely to be interested in the safest plant regard- 
less of cost. While these two extremes may provide useful information by bracketing 
the feasible design space, neither is likely to be considered a “good” design because of 
the conflicting nature of the objectives (i.e., an acceptable level of safety will not be 
achievable at low cost). Variations between these two extremes must be investigated, 
but “between” requires definition itself. Certainly the team is not interested in any 
design which is less safe and more costly than some other valid design. In fact, they 
are interested only in those for which no safer and less expensive valid design exists, 
i.e., those which are not dominated by any other feasible design. This, essentially, is 
the concept of Pareto optimality and is the only way that multi-criteria optimality 
may be defined without subjectivity.?? 

Although usually attributed to Pareto [69], this concept was proposed in an earlier 
work by Edgeworth [24]. The definition may be stated formally as follows [74]: 


“A feasible solution to a multi-criteria (multi-objective) optimization prob- 
lem is Pareto optimal if there exists no other feasible solution that will yield 
an improvement in the performance of one criterion without causing a de- 


crease in performance of at least one other criterion.” 


This means that there are multiple, and perhaps infinitely many, optimal solutions 
in the Pareto sense to most engineering design problems. These solutions are all 
equally “good”; without subjective preference information there is no distinguishing 
among them in terms of desirability. The set of all Pareto solutions is variously 
known as the Pareto set, the efficient set, the Pareto frontier, or the non-dominated 
frontier. It always lies at the boundary between feasible (realizable) objective space 
and infeasible (non-realizable) objective space, although the converse is not necessarily 
true, as shown in Figure 1-8. These plots represent hypothetical situations where 
every conceivable combination of design parameters has been evaluated (or built) and 
each combination’s relevant objective values (A and B) have been determined. The 
Mi Cost'and safety are obviously Rete only factors which would belcenseeredty Henideslenina aenucteanipomnn 


plant, but they are certainly primary factors. One difficulty in using the Pareto definition of optimality is how to 
interpret and display results when several objectives are present. 
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Figure 1-8: Notional Pareto frontiers in two-objective space, with objectives A and B to be mini- 
mized. (a) Continuous frontier, showing regions where objective B can improve significantly with 
little tradeoff of objective A. (b) Discontinuous frontier, with non-Pareto optima at the concave 
boundary of the feasible solution space. (c) Discontinuous frontier with two non-Pareto, feasibility- 
limited regions. (d) Non-conflicting objectives, showing local optima “trap” at the feasibility bound- 


ary. 


non-dominated frontier in each case makes up a portion of the feasibility boundary, 
obtained by plotting all of the feasible variations in objective (A-B) space. 

It is difficult to imagine how an algorithm could be designed to search for infinitely 
many (or even several) different solutions simultaneously. This is why optima-defining 
methods such as those previously listed apply some form of user-supplied subjectivity 
to reduce the number of candidates—in most cases to one. Depending on the method, 
this subjectivity will involve the assumption of one or more of the following: a global 
optimum and its location, the existence and values of independent optima for each 
objective, the mathematical formulation of a decision-maker’s preferences, objective 


thresholds and/or goals, or convexity of the Pareto set. These assumptions are ma- 
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nipulated such that the set of all objective values for any candidate solution can be 
combined into a scalar (for example, in the global criterion method, this scalar is the 
candidate's Euclidean distance in objective space from the location of the assumed 
global optimum). The combining of multiple objective values into a single measure 
of goodness is known as scalarization, and methods which combine objectives are 
known as scalarization methods. The value of the scalar provides the search process 
with a means of distinguishing among Pareto optima. The search through design (pa- 
rameter) space is thus guided toward a particular solution in objective space by the 
definition of the optimum. When and if this particular Pareto optimum is located 
by the search process, the corresponding design parameters are declared to be the 
solution. 

A serious objection to scalarization is that the form of the feasibility boundary 
is not taken into account. Accepting a definition of optima a priori amounts to an 
uninformed choice of one Pareto solution from among many. This implicitly assumes 
that gradients on the frontier are uniform, when in reality the frontier may contain 
concavities or even discontinuities. Rapidly changing gradients, or “knees,” in the 
frontier should certainly be considered in decision-making, as these are regions where 
a slight trade-off of one criterion may provide relatively large returns in others (the 
“regions of high return” in Figure 1-8a). This information is not available when using 
scalarization. Depending on the assumptions required for the chosen scalarization 
method, there is also a possibility that the search process will become trapped at a 
local optimum, such as that of Figure 1-8d.!2 None of the methods listed above can 
guarantee global Pareto optima [65]. 

Also, while the solutions defined by these methods are Pareto optimal if they ez- 
ist, they do not necessarily exist in an arbitrary solution space. For example, the 
assumption that a decision-maker’s preference function is a linear combination of 
objectives—a common feature among weighting methods—will not allow the identifi- 
cation of any optima in concave regions of the Pareto frontier, regardless of the search 
method used [65, 80].!° Even if the solutions are identifiable using the chosen defini- 
tion, the search mechanism may not be able to locate them. Most search processes 
rely on gradient information and therefore will have trouble with discontinuous or 
disjoint objective spaces, these being not uncommon in engineering design.'* Thus a 
multi-objective optimization algorithm can suffer from two general limitations: that 
" 12 Figure 1-8d also shows the only situation where subjectivity can be completely avoided in multi-objective opti- 
mization: compatible objectives with a single, feasible non-dominated solution. 

13 Note that concavity of the feasibility boundary does not necessarily indicate a discontinuous frontier, as may be 
seen by comparing Figure 1-8a and b. 


14 The gradient referred to here is the change in an objective value with respect to a design parameter, as opposed 
to the gradient of the frontier itself which is the change in one objective value with respect to another. 
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of its scalarization and that of its search mechanism. The degree to which these affect 
solution quality, of course, is dependent upon the nature of the particular objective 
space and cannot be generalized. 

Even if these limitations are acceptable, the user of a scalarization method must 
allow that any solution obtained is non-unique and is subject to change. Any new 
or altered assumptions, whether due to new information, a new decision-maker, or 
changing opinions of the same decision-maker will require that the problem be re- 
solved to locate a different non-dominated point. 

Based on these considerations, it seems that the ideal situation from a decision- 
maker or designer standpoint is to begin with a complete and accurate rendering of 
the Pareto frontier and proceed on the basis of this information. This has not gone 
unrecognized; there have been methods proposed which are capable of generating 
some or all Pareto solutions simultaneously for discrete problems (e.g., Rosenman 
and Gero {74]). Similarly, some authors of scalarization methods recommend that 
the assumptions be systematically varied and multiple solutions found prior to final 
selection (e.g., Dlesk and Liebman [20]). An algorithm which produces such informa- 
tion in a single pass, however, would seem preferable to one which requires repeated 
applications. Indeed, in light of all the issues discussed thus far in this section, the 
ideal multi-objective optimization algorithm would be capable of using the Pareto 
definition of optimality exclusively (thus avoiding a priori subjectivity), would lo- 
cate all Pareto optima simultaneously (thus exposing regions conducive to tradeoffs 
without repetition), and would not be hindered by pathological objective spaces. 

One approach which meets these criteria is a Pareto genetic algorithm. Genetic 
algorithms process populations rather than individual solutions and are therefore 
uniquely suited for defining a Pareto frontier in a single run. Since they do not 
rely on gradient information, they are unaffected by discontinuities and are not as 
likely to become trapped at local optima. This research attempts to exploit these 
characteristics in obtaining the frontier for ducted propeller submarines. The goal 
here will be validation of the method itself and demonstration of its superiority over 


scalarization, as well as drawing hydrodynamic conclusions from the results. 


1.4 Chapter Summaries 


Chapter 2 describes the submarine design and analysis program Ducted Propulsor 
Lifting Line (DPLL), which is eventually used as the evaluator for the genetic algo- 


rithm. The initial phase of this research involved an extensive re-write and upgrade 
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of DPLL, with the intention of making it suitable for use in full stern optimization. 

Chapter 3 provides validation of the revised DPLL by modeling published experi- 
ments and comparing the calculated results to measured data. The program’s thrust 
and torque calculations are verified, as well as the performance of a new boundary 
layer routine in predicting displacement thickness and flow separation. Sample out- 
put is included from a run typical of those in the optimization phase of the research, 
to provide qualitative verification of DPLL’s accuracy. Convergence characteristics 
of the program are presented. 

Chapter 4 includes a brief history of genetic algorithms and the theory behind 
them, along with a description of the algorithm developed during this research. The 
three alternative Pareto methods of selection used in this algorithm are described. 
Incremental and generational genetic algorithms are compared. The concept of com- 
petitive species in a genetic algorithm, apparently unique to this research, is intro- 
duced. | 

Chapter 5 documents the results and presents the Pareto frontier obtained by 
each of the three selection methods. The cumulative results of the GA are used to 
estimate the location of the dominated feasibility boundary. The composition of the 
entire feasibility boundary, in terms of design parameters, is presented. The results 
of a three-objective optimization are given, with minimal cavitation included as the 
third objective. 

Chapter 6 presents hydrodynamic and optimization-related conclusions. Several 
specific topics are recommended for future work, including possible improvements to 


DPLL and further submarine optimization issues. 
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Chapter 2 


DPLL version 2.0 


2.1 Description and Motivation 


Ducted Propulsor Lifting Line (DPLL) is an exploratory design and analysis code for 
underwater vehicles, in particular those with ducted propulsors. Using input geom- 
etry, design parameters, and operating conditions, it calculates all torques, thrusts 
and drags acting on the vehicle and its propulsion components (see Figure 2-1). Since 


these forces are assumed to be strongly coupled, the solution process is iterative. 


duct pressure drag 


rotor thrust <= duct viscous 
%A and induced drag 


duct lift 
body pressure drag Rowe 
stator viscous 


drag 


body viscous drag rotor viscous 
drag 





Figure 2-1: Coupled forces acting on an underwater vehicle. 


Version 2.0 of DPLL, developed during this research, is the result of an extensive 
re-write of Taylor’s original version [82] and includes several new features. Basic 
methodology remains largely the same as that described in Taylor’s written thesis [83]. 
Hereafter, version 2.0 will be referred to simply as DPLL. 

The distinguishing characteristics of DPLL are its ability to model an underwater 
vehicle as a system, thereby capturing the interaction among the various components 
involved, and to account for the effects of a contracting wake. It is also capable of 
converging to a self-propelled condition (i.e., balancing thrust against drag), predict- 


ing flow separation on the body, and estimating the effect of wake fraction. These 
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features make it more accurate and adaptable than simple parametric estimation 
such as [4] or non-coupled lifting line/actuator disk analysis such as [59, 93, 46] in 
the exploratory design phase. 

Like any model, DPLL is characterized and limited by its assumptions. In keep- 
ing with its intended use as an exploratory design tool, these assumptions are aimed 
at achieving the greatest accuracy possible while emphasizing simplicity and mini- 
mal run time. The use of a lifting line model for propeller blading, rather than a 
three-dimensional vortex lattice, is an example of a tradeoff between these objectives. 
Other simplifications include the use of assumed friction coefficients for calculating 
viscous drag and a zero thickness, mean camber line duct model. For these reasons, 
DPLL is less precise than state-of-the-art force/flow analysis tools (e.g., Reynolds- 
averaged Navier-Stokes, or RANS, routines coupled with propeller design codes such 
as PBD [56, 5, 6]).1. Such tools are computationally intensive and require consid- 
erable setup and run time (experienced users may require several days, given some 
arbitrary configuration, to obtain an accurate converged solution). In comparison, 
DPLL trades off high-order accuracy for greater speed, simplicity and ease of use. 

In short, DPLL is intended to fill a current gap between over-simplification and 
unnecessary precision in design at the exploratory level. It is motivated by the need 
for a relatively simple code with enough accuracy, flexibility and speed to be useful 


in preliminary submarine design. 


2.2 Methodology 


Propeller-induced wake vorticity in the DPLL model is concentrated into discrete 
streamtubes shed from lifting segment endpoints, as shown in Figure 2-2. The paths 
of these streamtubes (or streamlines, when represented in two dimensions) are deter- 
mined by the local velocity calculated at numerous intermediate points. These local 
velocities are affected by the presence of the body and duct as well as the vorticity 
contained in the streamlines themselves. Once their paths have been determined, the 
streamlines are discretized into free vortex segments. Propeller-induced velocities in 
the wake are then calculated by summing the effects of all bound and free vortex 
segments at numerous wake field points between the streamlines. The field points 
compose a propeller-induced velocity grid, which is interpolated or extrapolated to 
control point locations on the duct and body surfaces to give the propeller induction 


1The propeller blade design program PBD-14.3 and its earlier versions were developed under the supervision of 
Prof. Justin Kerwin at the MIT Marine Hydrodynamics Lab. 
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Figure 2-2: Centerline section of a notional DPLL wake, showing a typical grid point and the various 
components which affect velocity. No stator is present in this configuration. The rotor is modeled 
by five discrete lifting segments, each with its own input value of circulation. 


there.? Since normal velocity at control points is required to be zero—no flow may 
pass through the hull or duct—body and duct singularity strengths must be adjusted 
whenever new propeller-induced velocities are calculated. When such adjustments 
are made, the local, non-propeller-induced velocities in the wake are affected and 
streamline paths must be recalculated. This process is repeated until changes in grid 
velocities from one iteration to the next (among several other criteria) are within 
some tolerance. Convergence of the process results in an accurate total velocity field 


from which forces and torques on the lifting segments may be calculated. 


2.2.1 Lifting Line Model 


DPLL will model up to three propeller stages; these may be rotors or stators in any 
combination. The stages may have up to ten blades. Each blade is modeled by a 
lifting line, composed of vortex segments which induce rotational velocities normal 
to their axes in the surrounding flow. The magnitudes of these induced velocities are 
specified by the input circulation (T), or strength, of the segments. Figure 2-3 shows 
an example vortex segment oriented along the z-axis having length |r2 — x1|, zero 
thickness, and circulation l. At any point P in the vicinity, the velocities induced in 
the x and y directions (u and v respectively) are zero. The tangential component of 

2 Propeller induction is extrapolated rather than explicitly calculated at body and duct control points because 


these points lie on streamlines. Influence calculations become singular when the field point is on the axis of a vortex 
segment. This is also why propeller grid points in DPLL are located between streamlines. 
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Figure 2-3: Vortex segment parameters, from Kerwin [54]. 


induced velocity w is a function of the field point’s position relative to the segment 


and the circulation strength [54]: 


E 
wy) = Ae > a =] (2.1) 
l 
= 7 [cos 02 + cos 0] (2.2) 


Placed in a flow field, such a segment will not only induce velocities but will also be 


subject to a force according to the Kutta-Joukowski law: 
F=pUxT (2:8) 


In a lifting line propeller model, “blades” made of vortex segments are spaced sym- 
metrically and aligned radially about a common axis of rotation. In rotating through 
the surrounding velocity field and advancing forward with the vehicle, these segments 
experience forces in the plane of rotation and in the direction of the axis, creating 
torque and thrust respectively. 


3 As a consequence, 


DPLL assumes an ideal fluid and potential (zero curl) flow. 
vortex lines cannot begin or end anywhere in the flow field; any vorticity present on the 
lifting segments must be shed into the wake.* This allows calculation of streamline 
vorticity based on lifting segment strengths, which are DPLL inputs. Figure 2-4 
illustrates this relationship. A typical blade tip is shown, having a finite chord with 
an arbitrary distribution of bound circulation, yg. This bound circulation, essentially 
~ 3The only exception is when the hull boundary layer is calculated; see Section 2.3.1. 

1Kelvin's theorem of the conservation of circulation states that for an ideal fluid acted upon by conservative forces 


(e.g., gravity) the circulation is constant about any closed material contour moving with the fluid. Thus, any motion 
that starts from a state of rest at some initial time will remain irrotational for all subsequent times [67]. 
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Figure 2-4: Bound circulation related to free vorticity in wake for a general lifting surface and for 
the corresponding lifting line model, from Kerwin [54]. In practice, wake vorticity sheets tend to 
“roll up” into concentrated tip vortices. 


the velocity difference between the upper and lower faces of the foil, has dimension of 
circulation per unit length, or length per unit time. According to Kelvin’s theorem, 
the circulation around the isolated closed path in the upper part of the figure is zero. 
If the path is moved onto the blade and wake as shown without cutting any lines of 


vorticity, the total circulation remains zero and the following relationship must hold: 


ee / NS (2.4) 


Sle 


Pals) == / "atez)dos = =P (su) (2.5) 


Differentiating this with respect to s,, gives: 


dr, df 


ds,  dsw J, 


The same relationship between bound and free vorticity applies when the blade is 


7p(S2)ds2 = —7s(Sw) (2.6) 


w 


approximated by a lifting line. The chordwise distribution of circulation present on 


the actual blade is simply compressed to a vortex segment with the same value of 
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total circulation [(s,,). The two-dimensional blade surface, with finite chord and 
span, is replaced by a one-dimensional line (as shown in the lower part of the figure), 
with the assumption that the effect on field points off the blade surface is unchanged. 

If a lifting line is segmented and I is made constant (but not necessarily equal) 
on each segment, non-zero circulation gradients may occur at the nodes and the end- 
points. Lines of wake vorticity, or trailers, must begin at these locations and continue 
back through the wake to infinity in order to satisfy Kelvin’s theorem. The strengths 
of the trailers are known from Equation (2.6), and their paths are determined by local 
velocities in the wake as discussed at the beginning of this section. They will induce 
velocities on each other and on the bound lifting segments themselves, in accordance 
with Equation (2.1). The sum effect of all lifting segments and all trailing vortic- 
ity at a field point may be calculated numerically (by discretizing the free vorticity, 
as is done in DPLL) or analytically. Analytic methods require the assumption of a 
constant radius wake, however, and are therefore inappropriate for use in DPLL. 

A similar lifting line model is used for stators, which are non-rotating blades or 
guide vanes used to manipulate tangential velocities. Stators, like rotors, require input 
circulation distributions and are subject to torque and thrust. Stators may be placed 
upstream of the rotor to provide pre-swirl (with the intention of decreasing rotor 
torque), downstream of the rotor to recover wake energy in the form of tangential 


velocities, or both.? 


2.2.2 The H Function 


The circulation distribution of a lifting line propeller may be defined in terms of the 


velocity induced at the propeller plane by the sum effect of its helical trailers [46]: 


T(r) = a E (2.7) 


where U¿ is circumferential mean induced tangential velocity and Z is the number 





of blades. This quantity may be non-dimensionalized using a convenient reference 
length and velocity: 

P(r) 
Ven 


DPLL uses maximum hull circumference and forward speed as the reference quantities 





G(r) = (2.8) 


(the rationale for using body radius as opposed to the more traditional propeller radius 


5DPLL will also model multiple rotors, e.g., a contra-rotating configuration. This capability has not been tested 
or validated. 
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for normalization is discussed in [83]): 





Pr) 
= Z. 
CW) = a Re, 2.9) 
which becomes, using Equation (2.7), 
os 


The derivation of Equation (2.7) assumes an infinite number of blades uniformly 
distributed around the hub; its accuracy falls off rapidly for propellers of less than 
about ten blades [54]. If finite blade effects are to be taken into account, it is more 
appropriate to specify circulation in terms of local induced velocity (us) on the lifting 
lines, which is generally greater than the circumferential mean. This is the motivation 


for the H function proposed by Kerwin [55]: 





H Vat To 21 
== (2.11) 
Using the DPLL normalization quantities from above, this becomes: 
LA u” 
H(r) = 22 = G(r)— 2 
(r) = ER = air) (2.12) 


Lifting line segment strengths in DPLL are specified by H values at the hub and tip— 
the innermost and outermost segments, respectively—of each stage. These values are 


interpolated linearly to determine the circulation at intermediate points. 


2.2.3 Goldstein Factors 


Induced velocities in a propeller wake become circumferentially uniform downstream, 
regardless of the number of blades present. If local H values are used to specify 
bound vorticity strengths, as they are in DPLL, a relationship between local and 
circumferential mean velocity is needed in order to calculate the strength of wake 
vorticity using Equation (2.6). Such a relationship may be obtained by borrowing 
results from idealized propeller analysis. 

In 1929, Goldstein first solved for the self-induction of an optimum, finite-bladed 
lifting line propeller in uniform flow [38, 8]. The solution proved to be a function only 
of geometry—the number of blades, the pitch-to-diameter ratio, and the radius of shed 
vorticity. Its relationship to the infinite-bladed solution is known as the Goldstein 
reduction factor [54]: 


Ao) = ne (2.13) 


u(r) 
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where wu; is circumferential mean induced tangential velocity and uj is local induced 
tangential velocity. 

If a wake of constant radius is assumed, the Goldstein factor can be calculated 
directly using Equation (2.7) for a and the individual filament (single blade) equa- 
tions of Lerbs and Wrench [59, 93] for už. Unfortunately, the constant radius wake 
assumption is not appropriate in DPLL and neither of these calculations is suitable. 
However, it is reasonable to believe that the relationship between the two velocities is 
nearly the same regardless of wake contraction or the optimality of circulation. If that 
is so, then induced velocity at a lifting line can be approximately related to mean 
induced velocity in the plane of the propeller for contracting wakes. The velocity 
formulas mentioned above are simply used as if the configuration were ideal; the ratio 
of the two velocities is taken and used to convert known local induced velocity on a 
blade to unknown mean velocity in the propeller plane, or vice-versa.? 

Such use of idealized ratios to compute velocities in non-ideal configurations was 
proposed by Taylor [83]. He referred to these ratios as “generalized Goldstein factors” 
when describing their use in the original DPLL. The methodology remains the same 
in the current version. Non-dimensional local circulation values are used to determine 


the corresponding circumferential mean values via the generalized Goldstein factor k: 
lr nar) (2.14) 


where a value of x is calculated for each lifting segment once per iteration as if idealized 
conditions held. The G values thus obtained are used to determine the strengths of 


the trailing vorticity per Equation (2.6). 


2.2.4 Calculation of Propeller Induction 


The assumption of a circumferentially uniform wake simplifies the calculation of trailer 
paths and allows the wake to be represented in two dimensions, as in the centerline 
slice of Figure 2-2. The actual velocities, however, are induced in three dimensions. 
Each segment shown in Figure 2-2 represents a three-dimensional vortex ring, formed 
by sweeping the segment about the axis of symmetry. A typical, simplified ring is 
shown in Figure 2-5. The z and r coordinates of the two-dimensional segment’s 
endpoints are known from the growing and discretizing of the streamline. If the 
streamline were represented in three dimensions, however, it would be seen to follow 
a helical path; that is, the endpoints of the segment would be displaced tangentially. 
amet needed for the eenberanclWrenchiformulsaiwiichadoestnoteupeau inl & quaciona 2mpmel henna 


angle. This is the angle between the intersection of the chord line of the section and a plane normal to the propeller 
axis, assumed equal to the departure angle of the trailer at the lifting line. 
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(8) = constant 





Figure 2-5: Relationship between 2-D trailer segments and 3-D vortex rings. 


The tangential (@) coordinates are found by dividing the segment’s axial span (dz) 
by the latest calculated local axial velocity to obtain elapsed time. Local tangential 
velocity is then calculated and multiplied by this elapsed time to obtain the tangential 
span of the segment. A subroutine in DPLL written by Buchoux [9] uses the three- 
dimensional endpoint coordinates to convert the segment to its ring equivalent and 
then calculates the influence of the ring at any given field point. This process is 
repeated to obtain the influence of all trailing segments (rings) in the wake on all grid 


points; it comprises the bulk of DPLL’s processing time.’ 


2.2.5 Duct Modeling 


Propeller stages in DPLL are surrounded by a non-rotating duct of zero thickness, 
modeled by discrete vortex rings distributed along the mean camber line. Like the 
propeller lifting lines and the wake vortex rings, these singularities experience forces 
and induce velocities in the surrounding flow. Duct vortex rings are similar in effect 
and form to the free vortex rings used to model the wake (Figure 2-5), but have zero 
axial span. 

DPLL can treat the duct as a design or an analysis problem, according to a user- 
specified flag in the input file. In design mode, the total duct load and a preliminary 

T The sum effect of all trailing segments at some field point does not result in the total velocity there, as it does 


not include the effects of the body and duct singularities or the freestream. These additional effects are necessary for 
growing new streamlines, but are calculated separately. 
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camber line shape are input. The duct's overall angle of attack and scaling factors 
for the input shape are re-calculated at each iteration to result in the required load 
and zero normal velocities at duct control points. Providing the option of duct 
design complicates the program considerably; however, it also avoids some potential 
problems which will be discussed presently. 

Design mode requires the axial extents of the camber line as well as the load 
and an initial shape as inputs. The radial position of the duct is also specified, 
although indirectly, in that it must intersect the tip of a selected governing propeller 
whose radius is fixed. The camber line scaling factors, angle of attack and vortex ring 
strengths are computed given the input parameters and the latest velocity field around 
the duct. When the camber line is moved to its newly solved position, however, the 
mutual induction between the duct vortex rings and the body source rings is altered. 
The result is non-zero normal velocities at control points, meaning that the system 
must be re-solved, starting from the new duct position and using updated induced 
velocities. Duct design is therefore an iterative sub-process within DPLL’s main loop, 
involving repeated load computation and camber line adjustment until the calculated 
load is within some tolerance of the input. 

Figure 2-6 illustrates the freedom DPLL has in adjusting duct camber. A notional 
camber line which might result from an input B-spline vertex file (the preliminary 
camber shape mentioned above) is shown and denoted as “base camber.”? The lighter 
lines are a few of the shapes that this input camber line might be allowed to assume 
during the duct design process. These lines represent changes in camber magnitude 
and are produced by scaling the radial coordinates of the defining vertices, thus 
redefining the B-spline itself. 

Camber scaling is one of two degrees of freedom available to DPLL when re- 
designing the duct. The other is the duct angle of attack, involving a rigid rotation 
of the camber line about its pivot point (the pivot point is the previously mentioned 
intersection of the camber line with the governing propeller tip). For a given velocity 
field, base camber and design load, there is a unique combination of camber scaling 
and angle of attack which results in zero normal velocity at all duct control points.’ 

Alternatively, DPLL may be run in analysis mode. This simply causes the duct 
design routine to be bypassed; the camber line specified by the input B-spline vertices 
is never altered, except for an initial radial translation to intersect the governing 
 8In DPLL, total duct load is the algebraic sum of the non-dimensional strengths of all duct vortex rings. Duct 
control points lie on the camber line between rings. 

9 camber line with an inflection point is shown for purposes of illustration. In practical foils, the direction of 


curvature normally does not change. 
10This unique solution is not necessarily feasible; see Section 2.3.6. 
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Figure 2-6: Possible variations of an input duct shape due to scaling. The nodes shown are at the 
axial locations of the camber line’s B-spline vertices. The vertices themselves, to which the radial 
scalings apply, are not shown. 


propeller tip.'! The duct angle of attack is constant and implicit in the curve defined 
by the input vertices. Duct load is calculated rather than specified; it will fluctuate 
and settle, along with the surrounding velocity field, as the iterations progress. 

It is an apparent contradiction that an arbitrary camber line shape and angle of 
attack may be specified in analysis mode, while a unique combination of camber scal- 
ing and angle of attack exists for a given load in design mode. This is explained by a 
previously unmentioned constraint on DPLL’s duct design solution—the requirement 
of small leading edge load, or shock-free entry. The requirement of shock-free entry 
simply ensures that the leading edge of the camber line is reasonably well aligned 
with the local inflow. It involves the two most upstream vortex rings on the camber 
line. If the circulation values of these rings are [; and ['2, and the distances between 
the first two duct control points (i.e., the one between the rings and the one just 
downstream of the second ring) are Az, and Azz respectively, then the shock-free 


constraint is simply 


The vortex rings and control points are cosine spaced on the camber line, so the 
forward panel is always the smaller of the two. This makes I’; small and thus forces 
the leading edge toward the direction of the inflow. 


11 This initial translation occurs in design mode also. 
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In a real flow, a misaligned leading edge can cause a cavitation bubble or flow 
separation, depending on the degree of misalignment, the fluid velocity and the am- 
bient pressure. Currently, DPLL has no way of predicting cavitation or separation 
anywhere on the duct. The solution will therefore converge regardless of leading edge 
alignment, but will be invalid if the corresponding physical duct is subject to leading 
edge cavitation. The shock-free constraint on design mode simply helps prevent the 
user from specifying physically infeasible ducts. 

Clearance must exist between a physical, non-rotating duct and the blade tips. 
This clearance, known as the tip gap, is not explicitly modeled in DPLL; that is, 
circulation of practically any magnitude is allowed at the blade tips and the lifting 
lines intersect the duct camber line. A non-zero tip load is one of the advantages of 
a ducted propulsor, as discussed in Section 1.1. While DPLL will accept vanishing 
circulation at the blade tips (i.e., an input of zero H at the outermost lifting segment, 
as would be expected for a wide gap or an open propeller), this condition may cause 
negative propeller-induced velocities in the wake due to the grid extrapolations used. 
If the magnitudes of these negative velocities happen to be greater than that of the 
potential flow in the vicinity, negative total velocity results and the program will 
be unable to grow valid streamlines in the following iteration. DPLL is therefore 
best described as a small gap model, carrying load at the blade tips but assigning 


secondary importance to tip gap effects.?* 


2.2.6 Hull Geometry 


The hull profile is input as a set of B-spline vertices and modeled by submerged source 
rings. These rings are similar in concept to the vortex rings discussed above, but they 
induce radial rather than tangential velocities (a typical body source ring is shown 
in cross-section in Figure 2-2, along with the direction of its induced velocities). The 
strengths of these source rings are re-calculated during each wake iteration to zero 
the normal velocities at all hull and duct control points. The number and density of 
control points on the hull is controlled by the user; a source ring corresponds to each 


control point so that the system is properly constrained. 


2.2.7 Thrust and Torque Coefficients 


Thrust and torque coefficients in DPLL (Cr and Cg respectively) are normalized to 
reference velocity and body cross-sectional area for reasons discussed in [83]. This 


12 For an analysis of the effect of tip gap variation on propeller performance, see McHugh [63]. 
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is in contrast to the use of propeller rotation rate and diameter for normalization, 
resulting in the more common coefficients Ay and Kg. 

The total tangential velocity seen by a lifting segment is due to its rotation, its 
self-induction, and the induction of all other lifting segments present. For an incre- 
mental annulus of a Z-bladed lifting line propeller, the thrust produced according to 


Equation (2.3) is: 
6T = Zp(wr + ús)T ór (2.16) 


where ù; is the total local induced velocity and T is the circulation of the lifting seg- 
ment across the increment. Using the quantity 5pV? Ag as the reference force, where 
Ag is the cross-sectional area of the body at its maximum radius, and substituting 


the definition of G from Equation (2.9) gives the incremental thrust coefficient: 


tt oot 





(2.17) 


CT = — 
OG VG ( V = 


Thus the total thrust coefficient for the propeller is: 


Cr = 42 | - (A) mE (2.18) 


where Ry and Rr are the hub and tip radii, respectively. 

Likewise, the total axial velocity at a lifting segment is due to the combined effects 
of forward speed and induced velocities from the hull, duct and propeller stages 
(including self-induction). The torque on an incremental annulus of a lifting line 


propeller, by Equation (2.3), is: 
0Q = Zp(V, + ta)P rér (2.19) 


where ú, is the total local axial induced velocity. Using the quantity 3pV, ABD as 
reference torque, where Dg is the maximum body diameter, and again applying the 


definition of G, the incremental torque coefficient is: 


r or 





ü 
Cgo = 2Z | 1+ = ) G= 2.20 
va ( y za ORs Rp (2.20) 
and the total torque coefficient for the propeller is: 
Rr a 
Co = 2z | (: + “| G(r) dr (2-21) 
Pe Vs Rb 


where Ry and Rr are again the hub and tip radii, respectively. 
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2.2.8 Summary of Inputs 


Table 2.1 summarizes the primary inputs for the stand-alone version of DPLL, most 
of which have been discussed in the preceding sections. The modified version of DPLL 
used for the optimization portion of this research requires essentially the same inputs, 


but some are held constant. These will be noted in Chapter 4. 


DPLL Input 















]— 
propeller rotation rate JB ES (advance coefficient) 


circulation at blade root? oo T (non-dimensional local circulation)* 


H A (non-dimensional local circulation)* 


miscellaneous number of control points, lifting segments, and wake panels, 





















location of ultimate wake, friction coefficients, speed, etc. 


a One value required for each propeller stage present. 

b V, = forward speed, n = rotation rate, Rg = max body radius. |J| > 10.0 implies a stator. 
€ uł = local self-induced tangential velocity; Z is the number of blades. 

d N is the number of vortex rings used to model the duct camber line. 

e The governing propeller is an arbitrary designation when multiple stages are present. 


Table 2.1: DPLL input parameters and their format 


2.3 Major Revisions and New Features 


Following are the significant changes incorporated during the re-write of DPLL v1.0. 
New features include hull boundary layer modeling, hub vortex drag, self-propulsion, 
calculation of cavitation index, and several new convergence criteria. The duct design 


process is extensively revised and includes failure recovery and convergence criteria 
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of its own. General normalization procedures for both internal and external values 
are standardized throughout. Hereafter, version 2.0 will again be referred to simply 


AI 


2.3.1 Boundary Layer Modeling 


DPLL computes boundary layer displacement thickness, momentum thickness and 
shape factors at hull control points and uses these calculations to predict flow separa- 
tion and estimate the effect of wake fraction. It employs a boundary layer marching 
routine (hereafter, MRCHBL) originally written as a stand-alone program by Prof. 
Mark Drela!® and used by permission [22]. MRCHBL is run in direct mode, meaning 
that edge velocity—the freestream velocity just outside the viscous boundary layer, 
often given the symbol u.—is prescribed. This input comes from DPLL’s inviscid 
solution for velocities at the hull surface. 

It is an approximation to assume that the outer velocity in the viscous boundary 
layer is equal to inviscid surface velocities on the hull. The error involved is not 
prohibitive at DPLL’s level of accuracy, however, and this technique gives very good 
results to first order as will be seen in Chapter 3. In general, it is possible to iterate 
a process such as this, by returning the calculated boundary layer thickness to the 
inviscid solver, adjusting the source strengths appropriately, calculating new inviscid 
“edge” velocities at the displaced surface, and so on. This is known as a weakly or 
loosely coupled process.'4 The boundary layer routine in DPLL, however, is run only 
once, after all inviscid iterations are complete. This is because coupling will cause 
direct interacting methods such as MRCHBL to diverge if separation occurs [66], and 


separation must be expected at times when modeling full stern vehicles. 


Estimation of Wake Fraction 


Wake fraction is a beneficial phenomenon associated with propellers operating in 
viscous boundary layers. It results from the velocity reduction and resulting boundary 
layer growth near the surface due to friction and pressure effects. If the velocity 
parallel to the surface of the hub (i.e., normal to the lifting line) is reduced without 
affecting tangential velocity, then Equation (2.3) reveals that less torque is required to 
produce the same amount of thrust. This could be considered an efficiency “increase”; 
it is more precisely a partial recovery of energy lost to viscous friction. In any event, 
it is a condition which becomes more prominent as the fullness of the stern increases 
and which cannot be captured by the inviscid calculations of DPLL alone. 


13 Associate Professor of Aeronautics and Astronautics, MIT. 
14 A strongly coupled scheme solves the boundary layer and inviscid flow equations simultaneously. 
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Velocity within a viscous boundary layer increases from zero at the surface (the 
empirically determined “no-slip” condition [75]) to the freestream value at its edge, 


as shown schematically in Figure 2-7. A rotor operating within a boundary layer will 


hull surface 





Figure 2-7: Typical boundary layer velocity profile. 


therefore see a radius-dependent velocity reduction as compared to the corresponding 
inviscid flow. The velocity profile shape is in general dependent upon the local pres- 
sure gradient and the curvature of the surface. It is also a function of edge velocity, 
the viscosity of the fluid, and whether the flow in the boundary layer is laminar or 
turbulent. 

Laminar and turbulent boundary layer flows are governed by the same integral 
momentum equation, obtained by integrating a combination of the continuity and 


simplified z-momentum equations term by term across the boundary layer [21]: 


dé 0 due Tw _ Cy 


¡iS 5 (2.22) 
where 
§ = momentum thickness = =: [ — uz)uz dy (2.23) 
e Jo 
6“ = displacement thickness = — Sa — us) dy (2.24) 
e Jo 
H = momentum shape factor = = (2.25) 
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(2.26) 


Ty = Wall shear stress = py 2 
Oy 
— , 2Tw 
C; = skin friction coefficient = (2.27) 
p 


2 
Ue 


wW 





In laminar flow, the three variables 0, H and C's can be reasonably well related with 
one-parameter velocity profile approximations. For submarine-size vehicles moving 
at typical speeds, however, the transition from laminar to turbulent flow within the 
boundary layer occurs near the bow and the boundary layer remains turbulent over 
the remainder of the vehicle. Unfortunately, the shape of the velocity profile for a 
turbulent boundary layer is not as easily resolved. Many different correlations and 
empirical relations have been put forth to allow closure of Equation (2.22) as it applies 
to turbulent flow; there are in fact over 50 different proposals in the literature [90]. 

One of the simplest approaches is to assume a profile shape consistent with empir- 
ical data. A surprisingly accurate approximation, first proposed by Prandtl [72] and 
recommended for general use by White [90], is the 1/7-power law: 


Ue => (+) (223) 


where dg9, the boundary layer thickness, is the value of y at which there is a 1% differ- 
ence between the boundary layer velocity and the freestream velocity.** Substituting 
Equation (2.28) into Equation (2.24) and replacing the upper limit of infinity on the 
integral with 099 gives a necessary relationship for estimating the effect of turbulent 


wake fraction: 
099 — 86” (2.29) 


DPLL makes use of the power law to estimate the velocity profile of the boundary 
layer at the rotor. Upon convergence of the main duct/propeller/wake loop, invis- 
cid surface velocities at all hull control points are known. DPLL sends the control 
point coordinates, their inviscid surface velocities and the global Reynolds number to 
MRCHBL, which returns the corresponding values of 6* (among other parameters). 
The value of 6* at the rotor is extracted and used to find the displacement thickness 
there using Equation (2.29). Assumption of the 1/7-power law profile then allows 
estimation of velocity reduction at any radius (i.e., any y value) on the lifting line.** 
This reduction is calculated at each lifting segment control point, located midway 
Mele aiere o well-definednmeriacel beGeeen the boundary layer and the surrounding How, as the viscous effects 
theoretically extend to infinity. 


16The fact that the rotor lifting segments are always oriented normal to the local inflow and are therefore not 
necessarily exactly normal to the hull is not taken into account. 
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between the segment endpoints. Rotor torque, which during DPLL’s main loop is 
calculated using inviscid velocities, is then updated using the reduced viscous veloci- 
ties. Rotor thrust, which depends only on tangential velocities by Equation (2.3), is 
assumed to be unaffected by the hull boundary layer. 

An approach such as this is clearly a low-order approximation; its results must 
be evaluated with an understanding of its limitations and assumptions. For exam- 
ple, it must be noted that the power law profile is most appropriate for internal 
pipe flows and for flat plates with zero pressure gradient. The local radius of hull 
curvature in the vicinity of the propellers is reasonably large even for full stern sub- 
marines, which makes the flat plate assumption tolerable. The pressure gradient in 
the vicinity, however, is generally non-zero and in fact may be considerable due to the 
combined influence of the hull, duct and propeller(s). Favorable pressure gradients 
(dp/dx < 0) tend to “fatten” the profile relative to the power law approximation, but 
so slightly that the effect may be safely ignored at this level of accuracy.!” Adverse 
pressure gradients (dp/dx > 0), which are much more likely to be present near the 


stern of a submarine, have the opposite effect as shown in Figure 2-8. The normal- 
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Figure 2-8: Turbulent boundary layer data and the power law approximation. Huang data is from 
the 1976 experiments, afterbody 3, J = 1.07 and z/L = 0.915. 


17 See White [90] Section 6-4 for quantification of these and the following comments. 
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ized 1/7-power law profile is shown along with two sets of data. The first of these, 
from Huang [48], is taken on a propelled axisymmetric full stern body immediately 
upstream of separation. A strong adverse pressure gradient exists at the data loca- 
tion. The power law under-estimates the velocity reduction in the lower portion of 
the profile, but the error introduced by the approximation is not severe. The other 
data is from a 1954 experiment by Clauser [14], shown as the dashed line in the figure 
(the profile is actually a prediction from the theory of Mellor and Gibson [64], but 
represents the Clauser data extremely well). This data is taken in the presence of very 
strong adverse pressure gradient, presumably bordering on separation. The power law 
approximation is obviously quite poor across the entire boundary layer in this case. 
What is important to note from the figure, though, is that the change in the bound- 
ary layer profile is not proportional to the magnitude of the pressure gradient. The 
power law approximation, which is very good in the absence of a pressure gradient, 
is only slightly less accurate in a strong gradient such as that associated with a full 
stern vehicle (e.g., the Huang data). Only on the brink of separation does the model 
become inappropriate, as illustrated by the Clauser data. It is reasonable to conclude, 
then, that the use of this approximation in DPLL provides acceptable accuracy in all 
but the most extreme cases. For configurations which border on separation, some of 
which will be investigated during the optimization phase of this research, the effect 
of wake fraction will be under-estimated but should at least allow accurate relative 


comparisons to be made. 


Prediction of Flow Separation 


The boundary layer integral momentum equation (2.22) is obtained by combining the 
two thin shear layer (TSL) equations and integrating the result across the boundary 
layer to solve for skin friction. The TSL equations are (1) continuity and (2) a first- 
order approximation of the z-directed Navier-Stokes equation as the Reynolds number 
goes to infinity. These two equations can also be combined and integrated so as to 


solve for dissipation, resulting in a second differential equation [21]: 





dé* 0* du 
— = 2Cp - 3 2.30 
o o de g 
where 
0” = kinetic energy thickness = =| (u? — uz)uz di (2.31) 
Ue Jo 
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y dissipation = / T(y) = 
0 





Cp = dissipation coefficient = 


(2.33) 


€ 


If Equation (2.30) is divided through by a new parameter H*, defined as: 


x 


H™ = kinetic energy shape factor = > (2.34) 


and Equation (2.22) is subtracted from the result, the kinetic energy shape parameter 


relation is obtained: 

















A C} 0 du, 

A Den a Ta (239) 
dH dH H*[2Cp Cy 0 du. 
A et aes ee Va 


which is normally used with Equation (2.22) in lieu of (2.30) in two-equation bound- 
ary layer methods [21]. Given the inviscid velocities ue(ĉ) from DPLL’s converged 
solution, MRCHBL solves these equations for H and @ using internal correlations for 


ne Cs, and Cp: 


H = EN 
C's = Cs(H,0, y, Ue) 
o UE) 


If the calculated value of H exceeds a pre-specified turbulent limit—set to 2.7 in 
the DPLL implementation—separation is declared.!® If separation occurs upstream 
of the rotor, the inviscid flow field and force coefficients calculated by DPLL are 
to some degree invalid. This situation would justify rejection of DPLL’s solution 
when it is used in a stand-alone mode. In the optimization of DPLL described later, 
where relative rather than absolute criteria take priority, it is more efficient to simply 
penalize separated variants. This will be taken up in more detail in Chapter 4. 

Like the wake fraction estimate of the previous section, the prediction of flow 
separation by this method provides consistent results and allows accurate relative 
comparisons among variants. It does not match full viscous flow solvers (such as 
RANS routines) in terms of absolute accuracy, but runs several orders of magnitude 
faster and is adequate for exploratory design. Validation of the routine’s ability to 
predict separation and estimate wake fraction are included in Chapter 3. 

“ASIA. this event, MRCHBLcontimles withthe boundary oreraa onan ia a e 


values from points upstream and ignoring the input edge velocities. The invoking of inverse mode is taken to indicate 
separation in this research. 
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2.3.2 Hub Vortex Drag 


The drag due to the hub vortex contributes 
to total drag on the hull. A hub vortex occurs 
when circulation at a blade root is shed onto the 
hub in accordance with Equation (2.6). This 
vorticity generally follows a hull streamline for 
some distance downstream before detaching, in 


contrast to tip vortices which leave the blade 





surface immediately. Figure 2-9 shows a hub 

vortex produced by a stator acting alone. Figure 2-9: Hub vortex downstream of a 
Wang [87] determined that a Rankine vortex a 

structure is an appropriate model for a hub vortex. The velocity and circulation of a 


Rankine vortex of radius ro are given by: 


E ifr < 
a(r) = 4 m TS (2.37) 


Io otherwise 


2rr 


‘gr? ifr< 
A (2.38) 


Po otherwise 


and Coney [15] gives the drag force due to the vortex as: 


DE. (1 ate 3) r? (2.39) 
0 


16 rr 


where rp is the radius of the hub and To is the circulation on the blade at the hub. 
In DPLL’s hub vortex model, [9 is taken as that of the innermost lifting segment. 
When multiple stages are present, it is the algebraic sum of the innermost segments 
of all stages. The value of ry is empirical; a reasonable value is one-quarter the radius 
of the final (most downstream) propeller hub. 

Inviscid body drag in DPLL is determined by summing pressure forces over body 
panels. When the point of departure of the hub vortex is reached (where hull 
radius —+ ro), the summation is terminated and ambient pressure is applied to an 
imaginary panel normal to the z-axis at that point. The force due to the hub vortex 
from Equation (2.39) is then added to the truncated pressure summation to result in 


the total inviscid drag. 
2.3.3 Self-propulsion 


DPLL provides the option of floating the advance coefficient of the rotor, allowing 


the program to search for the rotor rotation rate which results in a self-propelled 
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condition—where the sum of the thrusts produced by all components of the propulsor 
balances all drags acting on the vehicle.*? If two consecutive iterations of DPLL's 
main loop result in an over- or under-propelled condition, the advance coefficient is 
adjusted accordingly and the iterations continue until total axial force on the vehicle 
is within some tolerance of zero.*° All other convergence criteria (see Section 2.3.5) 
remain in effect. Convergence is disallowed during the iteration immediately following 


an adjustment of advance coefficient in order to allow settling time. 


2.3.4 Cavitation Index 


DPLL performs a simple, parametric estimation of the rotor’s cavitation index. Cav- 
itation occurs when fluid pressure is reduced below its vapor pressure and therefore 
changes abruptly from a liquid to a gas (fluids are essentially incompressible, but 
cannot withstand any significant tension [67]). This is a fairly common phenomenon 
on the low pressure (suction) side of propeller blades and at blade tips; it is often ac- 
companied by increased propeller noise due to the eventual collapse of the cavitation 
bubbles. Cavitation can also result in erosion of blade surfaces [30]. The method used 
by DPLL to relate rotor operating parameters to a cavitation index is adapted from a 
correlation by Farell and Billet [25] and detailed in Appendix A. The cavitation index 


defined there is simply -C the negative of the minimum pressure coefficient in 


Pmin ? 
the flow. This minimum is assumed to occur at the blade tips, where the clearance 
flow between the tips and the duct interacts with the primary through-flow. Since 
the pressure coefficient for pressures lower than the ambient is negative, larger values 


of this cavitation index equate to increased likelihood and severity of cavitation. 


2.3.5 Convergence Criteria 


The main loop of the code is tested for convergence during each iteration. Simul- 
taneous satisfaction of all criteria terminates the iterations and sends the program 
to the boundary layer and post-processing routines. The convergence criteria are 


summarized below. 


Goldstein factors 


New generalized Goldstein factors are calculated for all lifting segments during each 
iteration. These values are then subjected to several constraints intended to prevent 

19 An alternative would be to float the circulation distribution on the blades, but this would be a much more 
complicated process when multiple stages are present. 


20 Consecutive iterations must both be either under- or over-propelled to trigger an adjustment. This prevents 
time-consuming oscillations about the self-propelled point. 
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unstable oscillations. The constraints place absolute as well as dynamic limits on 
the new factors, based on their values from the two previous iterations. If the newly 
calculated factors must be constrained, the constrained values are used for the re- 
mainder of the iteration but the initial unrestricted values are used for convergence 
checking. This convergence criterion is satisfied only when the unrestricted values 
of all segments are within some tolerance of the previous iteration’s values. For this 
research, the convergence criterion is a change of less than 2% in any Goldstein factor 


from one iteration to the next. 


Propeller induction on body 


Convergence of velocities in the wake is assumed to be directly related to the change 
in propeller induction at hull control points during consecutive iterations. These 
changes require adjustment of body source ring strengths so as to maintain zero 
normal velocity at the hull, and therefore indirectly affect the entire flow system. 
It is necessary that these induced velocities become relatively constant in order to 
declare convergence. For this research, the criterion for maximum change in propeller 
induction at any body control point from one iteration to the next is 1% of reference 


velocity (vehicle speed, V,). 


Duct load 


When DPLL is run in design mode, the duct load is analyzed at the beginning of each 
iteration by summing the non-dimensional strengths of all duct vortex rings. Since 
wake streamlines will have changed since the previous duct design, the analyzed 
load will generally differ from the design load even if the previous duct design was 
successful. The program proceeds with a new duct design, which modifies the camber 
line and angle of attack to re-acquire the input load. When the iterations have 
progressed to the point where the duct’s post-design load is within some tolerance of 
the pre-design analyzed load, this convergence criterion is satisfied. This essentially 
corresponds to stabilization of the velocity field in the vicinity of the duct. The 
tolerance used in this research is a change in non-dimensional circulation (Gauct) of 
less than 0.001. This is at most 0.25% of the input target for the range of loads used 


here. 


2.3.6 Duct Design Failure 


The duct design process adjusts the magnitude of the duct camber (not the basic form, 


which is input; see Section 2.2.5) and the angle of attack to result in the specified 


of 


load given the current local velocity field. The process may fail to converge for either 
of two reasons: large variation in the velocity field from one iteration to the next 
(causing the design process to produce a camber line which intersects the hull profile; 
this is more common in the early stages) or numerically valid but physically infeasible 


input parameters. 


Damped Convergence 


Failure of the duct design process is usually a result of the duct leading edge en- 
countering the body profile. If the duct design is having difficulty converging, DPLL 
will attempt to relax the design constraints until a valid solution is obtained. This 
is done by temporarily resetting the target (input) duct load to a value closer to its 
latest analyzed load. If the duct design then converges, the wake calculations proceed 
using the relaxed conditions and the original input load is attempted again during the 
next iteration, when the required camber line shift will presumably be less radical. 
Since the velocity field should settle as the iterations proceed, this allows “damped” 


convergence of some configurations which would otherwise fail. 


Invalid Configurations 


A valid solution does not necessarily exist for the duct design process. The unique 
solution for a particular duct load, propeller tip radius and velocity field may place 
the leading or trailing edge of the duct within the body envelope. This can be an 
unforeseen consequence of apparently valid input parameters. 

The version of DPLL used for the optimization phase of this research is slightly 
modified to handle this situation. It will attempt to increase the radius of the duct 
if the relaxation routine described above is unsuccessful. This involves a permanent 
change to the input governing propeller tip radius. Fortunately, this is not often 
necessary. When the situation does occur, it typically involves the combination of 
negative duct load, small propeller radius and a very full stern. Adjusting the duct 
radius is considered preferable to non-convergence (and certainly to run time errors) 
in these rare cases. The stand-alone version of DPLL v2.0 does not have this feature; 


it will simply fail when these conditions occur. 


2.3.7 Normalization 


All spatial values within DPLL are normalized to body length. The body is assumed 
to be axisymmetric on the horizontal z-axis and facing left; the bow and stern are 


at z/L = 0.0 and 1.0 respectively. A non-physical “sting” extends aft of the stern 
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through the ultimate wake to prevent infinite tangential velocities due to radial con- 
traction.** This is simply a small radius continuation of the body profile and is not 
considered part of the body length for normalization purposes. 

All velocities are normalized to the far field velocity, which is assumed to be purely 
axial and equal to the vehicle’s forward speed. Non-dimensional coefficients are gen- 
erally normalized using maximum body radius (or diameter, in the case of advance 


coefficient), which is itself normalized to body length. 





21The ultimate wake location is user-specified. Downstream of this location, the wake and its streamlines are 
assumed to be of constant radius and streamline discretization terminates. Induced velocities due to the portion 
of the wake extending from this location to +00 are calculated analytically, using the expressions of Hough and 
Ordway [46]. 
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Chapter 3 


Validation of DPLL 


This chapter documents the modeling of published experiments in DPLL and com- 
pares the resulting calculations with experimental data. Although the primary use of 
DPLL in this research will be for optimization—requiring accuracy only to the extent 
that relative comparisons can legitimately be made—the code is also intended for use 
as a stand-alone design tool. Some indication of its absolute accuracy is therefore 
needed as well. 

The validations to follow begin with DPLL’s calculation of thrust and torque, us- 
ing well-known hydrodynamic experiments. These are followed by testing of the pro- 
gram’s boundary layer routines, using published wind tunnel data. Sample flow field 
output and force calculations are then presented for a notional full stern submarine 
with a ducted multi-stage propeller, typical of those considered in the optimization 
phase of this research. This is included to qualitatively support the conclusion that 
DPLL results are reasonable and valid. Finally, simple convergence characteristics in 
terms of the number of body and duct control points are given. 

Unless otherwise noted, quantities in this chapter are expressed using the nomen- 
clature of the original experiments. This will usually require a conversion of DPLL’s 
normal output values. A complete list of symbol and variable definitions and descrip- 


tions may be found beginning on page 15. 
3.1 Thrust and Torque Calculations 


3.1.1 Ka-4-55 


A series of five screw propellers was used by van Manen at the Netherlands Ship Model 
Basin (NSMB) to investigate the effect of radial load distribution on the performance 
of ducted propellers [86]. These propellers were derived from the P/D = 1.018 
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variant of the K-4-55 series by varying the pitch distributions and profile sections.’ 
The results of this investigation led to the development of a new series having one 
of these five propellers (denoted as screw B in the report) as its reference. The new 
series became known as Ka-4-505. 

The five experimental propellers (A through E) were tested in combination with 
four different duct designs. Three of these designs, referred to as ducts 18, 19 and 20 
in the report, differed only in the curvature of their outer surfaces and only by a small 
amount. There was no reported difference among them in terms of performance. The 
fourth design, referred to as duct 19A, was of similar shape but had a flat outer surface 
and a relatively thick trailing edge. The performance of 19A was slightly worse than 
that of the other ducts. 

The report by van Manen documents the performance of these five propellers in 
duct 19 and the performance of the entire new Ka-4-55 series in both ducts 19 
and 19A.* From among these possibilities, the particular configuration of propeller B 


in duct 19 is chosen here for modeling (Figure 3-1). This choice is driven by three 





Figure 3-1: Ka-4-55 series reference propeller B with duct 19. 


l! Propeller series are designated by a group name, number of blades, and expanded area ratio (EAR; roughly the 
ratio of blade surface area to 7D? /4). The K-4-55 series, for example, is from the Kaplan (K) group, has four blades 
and an EAR of 0.55. 

2The report also presents cavitation as a function of thrust, and efficiency loss as a function of tip clearance for 
the original five propellers. Unfortunately, this level of detail is not captured by the DPLL model. 
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factors. First, the blade circulation distribution of propeller B, a necessary input 
for the DPLL model which is not documented for any of the propellers, has been 
calculated by McHugh [63] and is readily available. Second, van Manen's data is 
presented in graphical form and the resolution of the A-E data, given at a single 
P/D value in duct 19 only, is higher than that of the full Ka-4-55 series data, which 
is given at six P/D values in both ducts 19 and 19A. Third, the Ka-4-55 series data 
is primarily intended to compare the performance of duct 19 with that of duct 19A. 
DPLL’s simple camber line model cannot distinguish between these two ducts. 
DPLL normally operates on input values of circulation at the propeller hub and tip, 
interpolating linearly to determine values at intermediate lifting segments. However, 
a linear approximation will not model the data well in this case, as can be seen in 


Figure 3-2. DPLL’s source code is therefore temporarily modified so that circulation 
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Figure 3-2: Ka—4~55 propeller B blade circulation vs. non-dimensional radius, as calculated and 
modeled. The DPLL model uses six discrete lifting segments. 


at each lifting segment may be individually specified. Unfortunately, this modification 
requires the bypassing of Goldstein factor calculations (see Section 2.2.3). This is a 
limitation of the current version of the code. 

The experimental circulation distribution shown is not explicitly provided in the 


report. Rather, it is taken from an independent calculation by McHugh using a 
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coupled DTNS/PBD-14 analysis [63] and the given experimental blade, hub and duct 
geometry. The DTNS (David Taylor Navier-Stokes) algorithm, developed at the 
David Taylor Model Basin, is a member of the class of routines known as Reynolds- 
averaged Navier-Stokes (RANS) viscous flow solvers. As the name implies, these use 
the full incompressible Navier-Stokes equations to solve for a system's velocity field. 
In a coupled approach such as that used by McHugh, the velocity field obtained from 
a RANS solver is routed to a propeller code—in this case PBD— which computes the 
resulting forces on the blades.? The forces are then fed back to the RANS code, which 
re-computes the velocity field accounting for the new reaction forces in the fluid. The 
process is iterated until some convergence criterion is met. In this case, the converged 
solution provided the experimental blade circulation shown in Figure 3-2. 

The inflow velocity profile is not mentioned in the report, but is assumed here to 
be axial and uniform. To obtain this effect in DPLL, the hull of a long and slender 
pseudo-vehicle is used to represent the experimental hub and shaft. The lifting line 
blades and the camber line representing the duct are located at the axial midpoint of 


the hub-body, as shown in Figure 3-3. 


lifting line model 


mean camber line 





Figure 3-3: Graphical representation of Ka—4—55 as modeled in DPLL, showing use of long “body” 
to simulate the experimental hub and shaft and produce uniform inflow velocity at the rotor. 


Determing the correct duct parameters for modeling the experiment is somewhat 
involved, because DPLL requires that the camber line be defined by B-spline vertices. 
The use of B-splines is usually beneficial in that it consistently results in smooth 
curves; unfortunately, it is dificult to accurately reproduce a given curve without 
tedious trial and error. Also, minor changes in camber line vertex coordinates tend 
to have noticeable effects on duct thrust in the DPLL model. For these reasons, 


3The general method used by McHugh was documented earlier by Black [5, 6]. 


64 


attempts to input the absolute geometry of the Ka-4-55 duct with any degree of 
accuracy are unsuccessful, and DPLL is instead run in design mode (described in 
Section 2.2.5). Vertices producing a camber form visibly similar to that of the exper- 
iment are used, but more importantly the duct load is now input and varied so as to 
produce the experimentally measured duct thrust upon convergence. Some trial and 
error is required during this process, in order to obtain a converged solution which 
matches the experimental duct thrust, camber line and angle of attack. The results 


of the best match are shown in Table 3.1. The accuracy of the duct thrust coefficient 


Coefficient Experimental DPLL 


ae 0.1054 0.1007 
Kr total 0 041 1 0.3440 
Ko 0.0385 0.0370 


Table 3.1: Measured and calculated force coefficients for Ka-4-55 reference propeller, at J = 0.36. 


is irrelevant since it is intentionally matched. However, by matching experimental 
geometry and duct thrust, DPLL also calculates an overall thrust coefficient (due to 
duct and rotor) very near the measured value. Calculated rotor torque is within 5% 
of measured, but proves to be quite sensitive to small variations in the duct input 
parameters. Seemingly insignificant changes in the form of the input duct camber 
line can alter the calculated Kg by 5% or more. Calculated propeller thrust, on the 


other hand, is relatively unaffected by small variations in duct parameters. 


3.1.2 HIREP 


Zierke et al. used the High Reynolds Number Pump (HIREP) facility to investigate 
forces and flows in a multiple blade row machine [95]. This facility was installed in 
the test section of the Garfield Thomas water tunnel, located at Pennsylvania State 
University. A schematic of the experimental configuration is shown in Figure 3-4. 
The experiments were undertaken to provide a high Reynolds number database for 
the validation of numerical codes. This was motivated by the relative scarcity of such 
information, particularly for incompressible flows, and the numerous analysis codes 
under development. 

The design advance coefficient of the experimental setup was J = 2.31, as normal- 
ized to the rotor diameter. Tests were conducted at this and two bracketing values. 
Thrust and torque measurements were taken on the pump rotor; velocities and pres- 
sures were measured in the region of the pump rotor and stator. Measured velocities 


were the source for blade circulation distributions given in the report and used here. 
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Figure 3-4: HIREP experimental setup, showing the pump facility as installed in the Garfield Thomas 
water tunnel. 


The HIREP system is modeled in DPLL using a long, constant radius body to 
represent the hub in a similar manner to that used for the Ka-4-55 validation (see 
Figure 3-3, page 64). In this case, however, a horizontal duct camber line extending 
well forward and aft of the blades is used to represent the tunnel inner wall and DP LL 
is run in analysis mode in order to maintain its position. Only the pump stator and 
rotor are modeled in DPLL as propeller stages; the turbine stator and rotor, which 
served to turn the pump rotor in the experiments, are not considered here. 

Ás in the Ka-4-55 validation, blade circulation is forced to a known experimental 
distribution. The measured circulation on the pump rotor and stator blades is ap- 
proximated by six discrete values in DPLL, applied at the control points of the six 
segments making up each lifting line (see Figure 3-5). Note the attempt to average 
the measured circulation drop at the tip over the outermost lifting segment in the 
model. This is a crude way of approximating tip gap effects using DPLL’s small-gap 
model. 

A summary of the measured, design and calculated results at the design advance 


coefficient is given in Table 3.2. Measured and design torque were very similar; the 
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Figure 3-5: HIREP experimentally measured circulation (via laser Doppler velocimeter) and its 
approximation for input to DPLL. 


Coefficient Measured DPLL Design Point 
Kg 0.310 0.335 0.316 
Kr 0.890 0.806 0.780 


Table 3.2: Summary of measured, calculated and design force coefficients for HIREP at design ad- 
vance coefficient (J = 2.31). Experimental problems were encountered in measuring thrust. 


DPLL torque calculation is 8% greater than the measured value and 6% greater than 
the design value. Measured rotor thrust was 14% greater than the design value. This 
experimental discrepancy is not fully explained in the report, although it is mentioned 
that problems were encountered in performing thrust measurements [95]. A sampling 
window smaller than the period of shaft rotation and a pressure differential across the 
hub, which may not have been included in the measurements, are cited as possible 
causes. In any case, the thrust calculated by DPLL is 3% greater than the HIREP 
design thrust and within 10% of measured thrust, falling between the two values. If 
the reported measurement problems did in fact result in inaccurate readings (as the 
text of the report seems to indicate), then design thrust must be taken as the more 
accurate. The small difference between measured and design torque lends support to 


this assumption. In that case, DPLL’s torque and thrust calculations are both quite 
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close to the actual experimental data. 


3.2 Boundary Layer Modeling 
3.2.1 Shear Stress and Separation 


Huang et al. performed wind tunnel tests on three axisymmetric bodies having iden- 
tical forward offsets but different stern geometries [48]. The objective of these tests 
was to obtain boundary layer characteristics and pressure distributions for compari- 
son with theoretical calculations. The data obtained have since been widely used for 
validation of hydrodynamic and aerodynamic analysis codes. An example of the ex- 
perimental models is shown in Figure 3-6 (note that in this and the remaining figures 
of this chapter, the axial coordinate is x rather than z in accordance with the Huang 


reports). The models utilized a common nose cone but different mid-body and stern 


L/D = 10.975 





Figure 3-6: An example of the slender bodies used by Huang. Taken from Huang et al., 1976 [48]. 


sections. The stern sections, shown overlaid in Figure 3-7, varied in their rate of taper 
and were all of different lengths. Overall model length (3.066 m) was kept constant 
by changing the length of the parallel mid-body section as necessary. The tests rel- 
evant to the current research were performed at a Reynolds number of 5.9 x 10°, or 
air speeds of 28.9 m/s, assuming air viscosity of 1.5 x 107% m*/s. Therefore the data 
and calculations to follow should correspond, by a similitude argument, to a full-size 
submarine at a bare crawl or a 3 m autonomous underwater vehicle (AUV) at about 
2 m/s, or 4 knots. 

The stern section designated as afterbody 1 was slightly tapered and designed for 
a smooth transition of flow to the wake. Afterbody 2 was shorter with more abrupt 
taper and was intended to verge on separation under the experimental conditions. 
Both afterbodies 1 and 2 had profiles defined by polynomials. Afterbody 3 was the 
shortest and most abrupt of the three. Unlike the others, its profile was defined by 


a cosine curve; it included an inflection point aft of the shoulder and was intended 
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Figure 3-7: Afterbody geometries from Huang experiments (1976). Note the axis scaling. 


to cause flow separation. All three afterbodies were adapted to accept a common 
propeller hub, which began at x/L = 0.97. Shear stresses were measured on the body 
surfaces with and without a propeller installed and operating.* 

This experimental setup is modeled in DPLL by developing B-spline vertex files 
to result in the three profiles and supplying the correct Reynolds number for the 
boundary layer calculations. For the propelled body tests, the propeller geometry 
and circulation given in the report are input also. Again, the circulation distributions 
are forced and Goldstein calculations are therefore bypassed. 

These experiments did not involve ducted propellers. Since DPLL requires a duct 
if a propeller is present (propeller tip locations in the program are specified as a 
fraction of duct chord), the propelled tests are approximated by using a very short, 
zero load duct. This is intended to represent a non-ducted configuration as accurately 
as possible. 

In all Huang tests, a tripwire was used to induce the transition from laminar to 
turbulent flow at 2/L = 0.05. This effect is easily modeled, as MRCHBL allows the 
input of a manual trip location. As the routine begins marching downstream from 


4The propeller was driven by a 9-kW high-speed motor mounted inside the stern sections. 
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the bow of the vehicle, it continually checks for both the location of the manual trip 
and the first occurrence of conditions which would induce a natural transition to 
turbulence. If either is encountered, turbulent equations are substituted for laminar 
and are used thereafter. Thus an input trip location takes precedence if it is located 
prior to the point where natural transition would occur; otherwise, the trip location is 
irrelevant. Some test runs with the DPLL model indicate that the natural transition 
location for the Huang experiments would have been just downstream of the trip. 
Thus the manual trip takes precedence and natural transition plays no part in the 
experimental data or the calculations documented here. 

Shear stress data taken from the Huang report is compared to the corresponding 
DPLL calculations in Figures 3-8, 3-9 and 3-10. For afterbodies 1 and 2, differences 
in shear distribution between the propelled and unpropelled conditions are noticeable 
in both the data and the calculations (the increase of shear stress on the hull due 
to propeller operation is a well-documented phenomenon known as the thrust deduc- 
tion). Note that for both of these afterbodies, the propeller effect extends upstream to 
approximately r/L = 0.9, a distance of more than three propeller radii (the propeller 
has a radius of r/L = 0.03 and is located at z/Z = 0.983 on all three bodies). This 
observation supports the assumption that propeller suction can allow some control 
over the boundary layer. Also note that this effect is due to an open propeller; a 
ducted propeller might be capable of producing a more pronounced effect. 

Huang reported that no measurable difference existed between the propelled and 
unpropelled stresses on afterbody 3. As Figure 3-10 shows, the experiment resulted 
in flow separation at x/L = 0.92 (separation is indicated by shear stress dropping to 
zero). The shear distribution calculated by MRCHBL matches the data quite well, 
and inverse mode is invoked (thus predicting separation in the model) at /L = 0.93.° 
Since boundary layer calculations are done only at body control point locations, 
which in this case are spaced at intervals of approximately L/100, MRCHBL actually 
predicts separation between 0.92 and 0.93. 

The stress spikes in the calculations for afterbodies 2 and 3 at z/L = 0.84 and 0.87 
respectively are in agreement with experimental observations. On page 35 of Huang’s 
report, Figure 8 shows measured pressure distributions on the three afterbodies [48]. 
In all three cases, the form of the measured pressure distribution is very similar to the 
form of the stress calculations by DPLL, including the spikes on afterbodies 2 and 3. 
Unfortunately, these features are not reflected in the data due to low resolution. The 
experimental shear distributions of Figures 3-9 and 3-10 taken alone are therefore 


5Shape factor, rather than shear stress, is the criterion used in this research to predict separation but the two 
indicators are strongly correlated (see Section 2.3.1). 
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Figure 3-8: Measured and calculated shear stress on Huang afterbody 1 (1976) 
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Figure 3-9: Measured and calculated shear stress on Huang afterbody 2 (1976). 
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Figure 3-10: Measured and calculated shear stress on Huang afterbody 3 (1976). 


somewhat misleading, in that they do not capture the full magnitude of the pressure 
rise. 

Overall, the correlation between measured and calculated boundary layer behavior 
is fair for all three bodies in terms of absolute accuracy. The range of propeller 
influence is captured quite well, judging by the location at which the propelled stress 
diverges from the unpropelled stress in the data and in the DPLL calculations. Also, 
the relative differences between propelled and unpropelled stresses as measured and 
calculated are quite similar. Of particular importance to this research is the fact 
that DPLL correctly predicts the occurrence of flow separation on afterbody 3. The 
calculated location of separation also appears to be reasonably accurate. 

The experiments on afterbody 3 present a situation very relevant to this research, 
namely a configuration which results in separated flow under the influence of an open 
propeller. Of particular interest is whether surrounding the propeller with a duct can 
affect separation in the DPLL model. To investigate, a duct with significant load and 
arbitrary camber is added to the model of the propelled afterbody 3. All other inputs 
remain unchanged. The final geometry of the converged solution (using design mode) 


is shown in Figure 3-11. The resulting shear distribution is shown in Figure 3-12, 


T2 


overlaying the non-ducted calculations and data. Shear on the notional ducted system 
still decreases rapidly near the stern but no longer reaches zero; separation appears to 
be averted by the addition of the duct. Shape factor calculations confirm the absence 
of separation; that is, the shape factor does not exceed the prescribed turbulent limit 
of 2.7 at any location forward of the rotor. These results support the theory that 


ducts may be useful in avoiding separation on full stern submarines. 


3.2.2 Displacement Thickness 


Huang et al. later performed another set of wind tunnel tests on afterbodies 1 and 
2 [47]. The purpose was to measure pressure and velocity distributions across the 
boundary layer, allowing calculation of displacement and momentum thicknesses. 
These tests, like the previous ones, were performed at a Reynolds number of approx- 
imately 5.9 x 10°. Afterbody 3, which was known to cause separation at this speed, 
was not tested. 

Turbulence was again artificially induced at z/Z = 0.05. A complicated argument 
in the report seems to indicate that this trip location should actually be modeled 
as having occurred at z/L = 0.015. This latter value is used in the DPLL model; 
the difference this location makes in the results is insignificant. Figures 3-13 and 
3-14 compare the displacement thicknesses (6*) measured by Huang to the results of 
MRCHBL for afterbodies 1 and 2 respectively. The calculations are quite accurate 
in both cases, the worst discrepancy being a slight over-estimate very near the stern 
of afterbody 2. The near-separated flow conditions on this model may be partially 
responsible for the inaccuracy; MRCHBL solutions generally become very sensitive 


near separation. 


3.3 Sample Output 


The experiments modeled in the previous sections, while useful for purposes of val- 
idation, provide little qualitative information regarding DPLL’s performance in its 
intended role of submarine design. The HIREP and Ka-4-55 experiments, for ex- 
ample, involve isolated propellers and lack the effect of hull interaction. Both of the 
Huang experiments are at a very low Reynolds number relative to typical submarines 
and the propellers are not compound or ducted. 

This section is intended to lend qualitative support to DPLL’s validity, by doc- 
umenting the analysis of a notional full stern submarine with a ducted, multi-stage 


propeller. The inputs used for this run are shown in Table 3.3. The duct is slightly 
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Figure 3-11: Duct and wake geometry of notional propulsion system, from DPLL output files. 
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Figure 3-12: Shear stress on Huang afterbody 3 with notional ducted propeller. 









afterbody | 
— — — — calculated & (MRCHBL) 
= measured 5 (Huang) 


0.06 





0.04 


rL 


0.02 


0.7 0.8 0.9 l 
x/L 


Figure 3-13: Measured and calculated displacement thickness on Huang afterbody 1 (1979). 
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Figure 3-14: Measured and calculated displacement thickness on Huang afterbody 2 (1979). 
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Input Parameters 


duct load (G) rotor tip radius (%L) 0.020 
duct/rotor intersect (% duct chord) duct chord (%L) 0.069 
rotor tip circulation (H) stator tip circulation (H) —0.035 


rotor hub circulation (H) stator hub circulation (H) —0.031 
Reynolds number stator/rotor blades 5/5 
ultimate wake location (%L) lifting segments per line 5 





Table 3.3: Summary of DPLL inputs for notional submarine. 


decelerative, meaning that the axial component of its lift should contribute to drag 
rather than thrust. A stator is employed forward of the rotor to provide pre-swizl; 
like the duct, it should contribute to viscous and non-viscous drag. The specified 
Reynolds number for this run is 1.0 x 10°, which corresponds to an 80 m submarine 
traveling at 12.5 m/s, or about 24 knots.® 


Calculated force and torque coefficients are listed in Table 3.4 (recall that forces in 


Force Coefficients Propeller Coefficients 


body pressure force 0.036 | stator torque —0.049 
body viscous force  —0.150 | stator thrust —0.040 
duct axial lift —0.045 | rotor torque 0.035 
duct viscous force —0.013 | rotor thrust 0.221 
hub vortex drag —0.002 | final rotor speed (Jp) 1.75 
total drag —0.172 
total thrust 0.181 


total force 0.007 


Table 3.4: Force and propeller coefficients for notional submarine. 


























DPLL are normalized by SOV¿Ap, torques by 5pV¿ABDp). A positive force coeffi- 
cient indicates a force on the vehicle in the —z direction; a positive torque coefficient 
indicates a force on a component of the vehicle in the +0 direction. Body pressure 
force is small, but non-zero due to the induction of the duct and propeller on the 
stern. The duct contributes slightly to non-viscous drag, as expected due to its small 
negative load. Likewise, the stator produces a small negative thrust as a consequence 
of imparting pre-swirl to the rotor inflow. The force due to the hub vortex is negligi- 
ble, as much of the vorticity shed onto the hub by the stator is removed by the rotor. 
The net force on the vehicle is not precisely zero; this is simply due to the tolerance 
used for satisfaction of the self-propelled criterion. The final advance coefficient, hav- 
ing been automatically adjusted to achieve self-propulsion, is 1.75 as normalized to 
maximum body diameter. This corresponds to a rotor rotation rate of about 130 rpm, 


6 The configuration modeled in this section is one of the frontier variants from Figure 5-2. It does not correspond 
to any particular full-size or model submarine. 


given the forward speed of 12.5 m/s and assuming a full-size body diameter of 10 m. 
Figure 3-15 shows the tangential and axial velocity due to propeller induction 


alone. The flow region affected by the propeller is bounded above by the duct stream- 
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Figure 3-15: Propeller-induced tangential velocity (top) and axial velocity (bottom) for notional 
submarine. 


line and below by the hull streamline. The locations of the lifting lines representing 
the stator and rotor are clearly visible in the upper tangential velocity plot. The 
pre-swirl induced by the stator is quite apparent; it increases as the hull contracts 
between the stator and the rotor. Not all of this velocity is removed by the rotor. 
This is not necessarily non-optimal in general, as it is possible that viscous losses 
incurred in removing all tangential velocity from the wake are greater that the real- 
izable gains [88]. In the lower plot, only the location of the thrust-producing rotor 
is apparent since the stator induces no axial velocity on the flow.’ It is interesting 
to note that the suction effect of the rotor extends upstream a distance of approx- 
imately two propeller radii, a result which agrees with the data and calculations of 
Section 3.2.1. 

Figure 3-16 shows a schematic of the configuration along with total axial and 
radial velocity in the vicinity of the propulsor, where the effects of the body, duct and 


7The z-component of the total velocity vector certainly changes when the flow encounters the stator, but the 
prop-induced component of axial velocity does not. 


fe 


inflow (forward speed) have been included. The decrease in axial velocity forward of 
the rotor in this plot is quite pronounced, especially near the hull, as the total velocity 
vector becomes more radial there due to hull contraction. The location of the duct 
is clearly visible in the axial velocity plot, and it is obvious from the velocity jump 
across the camber line that the duct is subject to a force directed away from the hull. 
The axial component of this force, from the orientation of the camber line, appears 
to be directed toward the right (the +z direction, thus contributing to drag). This 
is verified by the calculated duct force coefficient. The radial velocity plot shows the 
effect of the duct’s slight negative load in “squeezing” the flow between the leading 
edge and the hull. There is little change in radial velocity across the camber line, as 


normal velocity across the camber line is zero and the angle of attack is small. 


3.4 Convergence Properties 


The full stern submarine of the preceding section is also tested for convergence by 
varying the resolution of the input geometry. During convergence testing, a self- 
propelled condition is not pursued—this prevents inconsistencies in final values due 
to tolerance in the force balance criterion. Instead, the program is simply run until 
calculated values are invariant from one iteration to the next. The fact that this 
results in a non-zero force balance for all tests is not relevant; the goal is to determine 
the minimum resolution at which calculations are reliable. 

Calculations appear to be insensitive to wake discretization. DPLL allows specifi- 
cation of four wake discretization parameters. Three are axial—the number of wake 
stations forward of the duct leading edge, the number between propeller stages, and 
the number downstream of the duct trailing edge. The fourth is radial—the number 
of lifting line segments (and thus the number of streamlines in the wake). None of 
these variables appear to have a significant effect on the solution, even when they 
approach the minimum values allowed by DPLL’s input-checking routines.® 

Likewise, the solution appears to be unaffected by the location of the ultimate 
wake as long as it remains downstream of a critical value. Advancing the location 
of the ultimate wake upstream from a baseline location of z/L = 1.10 produces very 
little change until the point z/Z = 1.02 is reached. Any location closer to the stern 
than this causes the program to fail. That such a close location allows accurate 
results is not entirely surprising; as seen in Figure 3-16 the wake streamlines of this 


8 Valid results are obtained for this variant with as few as 5 wake stations forward, 5 between and 10 downstream 
of the propellers. The number of lifting line segments has very little effect at values greater than 5, which is the 
minimum number allowed by DPLL. 
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Figure 3-16: Total flow field for notional submarine. The upper schematic shows the locations of the 
hull surface, duct camber line, lifting lines for the stator and rotor and their wake streamlines. The 
middle plot shows total axial velocity in the same region; the lower plot shows the corresponding 
total radial velocity. 
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configuration attain constant radii very soon after clearing the duct trailing edge. 
Discretization of body and duct control points—and therefore the source and vor- 

tex rings representing them—does have a significant effect on the solution. Figure 3-17 

shows the body pressure force and the rotor thrust as functions of body and duct res- 


olution. The total number of control points on the body, like the stations in the 
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Figure 3-17: Convergence of body pressure force and rotor thrust, as functions of body control points 
and duct vortex rings. Rotor thrust is calculated to three digit precision. 


wake, must be specified by three values—the number forward of the duct, between 
propeller stages, and downstream of the duct. In these tests, each of these values is 
increased in proportion with the total number of points on the body. Fluctuation at 
low resolution and good consistency as resolution increases is apparent in both plots. 
Other calculated values, in particular the torque on both the stator and the rotor, 


are also affected by these variations but to a lesser degree and are not shown. 


3.5 Summary 


Based on the results documented in this chapter, DPLL’s accuracy in force and torque 


calculations and in the modeling of the hull boundary layer is considered sufficient to 
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allow its use as an evaluator for optimization. T'hrust calculations are quite consistent 
with experimental measurements, particularly when the simplifications of the model 
are taken into account. Torque calculations, while somewhat less consistent than 
those for thrust, are near the experimental data in all cases and can be expected to 
give accurate relative comparisons among notional configurations. The prediction of 
flow separation, a crucial aspect of exploratory full stern design, appears reliable as 


does the calculation of boundary layer displacement thickness. 
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Chapter 4 


Pareto Genetic Algorithms 


4.1 Overview and Motivation 


As mentioned at the end of Section 1.3, a Pareto genetic algorithm is the search mech- 
anism chosen here for optimization of the DPLL model. Genetic algorithms (GAs) 
are stochastic, non-linear optimization routines loosely based on theories of biological 
evolution. In contrast to more traditional optimization methods which use gradient 
information to move between successively better points in solution space, GAs oper- 
ate on populations of solutions using models of natural selection, reproduction and 
mutation. As a consequence, they are uniquely adaptable to multi-solution problems 
such as the defining of Pareto frontiers. 

The goal of a GA is to continually evolve its population in the direction of im- 
provement. In single-objective and scalarized multi-objective GAs, improvement is 
measured by either the population’s average fitness and/or its most fit member; in the 
less common Pareto GAs, by the number of non-dominated solutions and their distri- 
bution in objective space.’ In all types, natural selection, or “survival of the fittest,” 
is simulated by giving preference in breeding to those members of the population who 
more closely resemble the assumed optima. Offspring are produced by mixing the 
design parameters (genes) of these preferred members, with some small probability 
of mutation during the reproductive process. The offspring are then evaluated for 
fitness and inserted back into the population. These steps are repeated until a ter- 
mination criterion is met, such as some number of generations without improvement 
or attainment of a pre-defined objective goal. 

Although rigorous proofs of their convergence characteristics remain somewhat 
elusive and controversial, GAs are increasingly popular in academia and have begun 
to see application in industry. They are inherently attractive, perhaps due to their 


1 Fitness in a genetic algorithm (or in the biological sense, for that matter) may be thought of as the relative ability 
to survive and produce offspring. 
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simplicity and novelty, and are therefore prone to indiscriminate usage. It is advisable, 
when faced with a practical optimization problem, to ensure that a GA is not chosen 
simply on the basis of its subjective appeal. This requires a basic understanding of 
what GAs are and are not, and of what they can and cannot do.? 

GAs are not reasonable alternatives to deterministic methods of optimization. Lin- 
ear and non-linear programming techniques, if applicable, will out-perform the best 
GA. Neither are GAs necessarily more precise or efficient than other non-deterministic 
methods, such as simulated annealing. In fact, their performance in certain test cases 
has been shown inferior even to simple stochastic hill-climbing [52]. They are not well 
suited for applications that require guaranteed response times (such as such as real- 
time control systems); their response time variance is in fact quite high in relation to 
other stochastic methods [62].° 

Rather, the distinguishing characteristics of GAs are their adaptability and ro- 
bustness; they are not hindered by pathological objective (cost) functions.* They 
are capable of optimizing discontinuous, disjoint, and other poorly-behaved functions 
and are one of the few alternatives available when no information is available a priori 
regarding the form of the solution space—in experimental optimization, for example. 

In searching an unexplored solution space, a tradeoff must be established between 
two conflicting objectives: exploiting the information contained in the best known 
solutions and exploring new regions. Pure exploitation makes exclusive use of exist- 
ing information; pure exploration abandons all known solutions on the premise that 
better solutions exist elsewhere. Hill climbing, or perturbation, routines are good 
examples of exploitative methods where the best known solution is always used as a 
starting point for improvement. Such techniques take a very narrow view of where 
opportunities for improvement lie; they are extremely dependent upon their starting 
locations and are susceptible to becoming trapped at local optima. At the other 
extreme, totally random search processes provide excellent exploration but make no 
use whatsoever of acquired knowledge. A properly implemented GA strikes a balance 
between these two extremes by simultaneously processing information obtained from 
many potential solutions dispersed throughout the solution space. 

The use of populations provides GAs with the ability to locate multiple optima 
—Z Actually, GAs-are capable of optimizing nearly andhing, with the exception of anomalous (decent emia: one 
which can confound their mechanisms. The question is really how well they perform on a given problem relative to 
other methods. Goldberg [32] provides an introduction to GA deceptive functions. 

3They have been shown, however, to be quite suitable for control system design [13, 12]. 
4The GA literature often refers to the evaluator as the “objective” or “cost” function, although it often is not a 


function in the mathematical sense and usually has nothing to do with monetary cost. Evaluation is discussed in 


Section 4.3.2. 
5 These examples and this method of describing a genetic algorithm are taken from Booker [7]. 
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simultaneously in the presence of single or scalarized multiple objectives. It also 
makes them excellent candidates for defining Pareto frontiers, alleviating the need for 
numerous independent runs to obtain the desired frontier resolution. 

Despite being burdened at times with overly optimistic expectations, GAs have 
demonstrated excellent utility in practice. Recent successes have been documented 
in such diverse applications as the tuning of power system stabilizers, scheduling 
of road projects and control of spacecraft rendezvous [71, 73, 53]. Multi-objective 
applications are now beginning to appear in the literature; examples include the 
design of broadband microwave absorbers, scheduling the production of chilled ready 
meals and structural synthesis of VLSI circuits [89, 81, 1].* 

Due to the stochastic nature of GAs and their relatively recent introduction, there 
are few absolutes in the implementation of basic GA operators. This is particularly 
true in the case of Pareto GAs, for which no documented database of results has been 
established. Implementation decisions during the design of a GA may include, for 
example, the sizing of the population, the proper mutation rate and the methods by 
which individuals are chosen to breed, live and die. General implementation guidelines 
based on published experimental results are scarce and at times contradictory even 
for non-Pareto GAs; scarcer still are analytical results and recommendations. This 
research is intended to supplement the Pareto GA performance data by comparing 
Pareto implementations. In particular, the effect of varying the selection method is 
investigated, as this is considered the most critical to performance. Variations will be 
compared based on their ability to locate and define the Pareto frontier for full stern 
submarines, using DPLL as the evaluator of propulsive efficiency and hull volume. 

The Pareto GA developed and tested here operates on populations consisting of 
multiple, competitive species. The need for a multiple species population is due to 
the different propulsor configurations allowed by DPLL; as will be seen, these have 
non-compatible characteristics and cannot be permitted to interbreed. The multi- 
species concept is apparently unique to this research (a multi-sexual, multi-objective 
GA has been proposed, but involves no competitive aspect [60]). Further discussion 


of the concept and its implementation will be taken up in Section 4.4. 


4.2 The Schema Theorem 


The beginnings of numerical optimization using models of biological processes can be 
traced to the late 1950’s, although credit is usually given to John Holland for pioneer- 


© Documented multi-objective applications generally employ scalarization. The microwave absorber application by 
Weile, Michielssen and Goldberg is one of the few Pareto implementations [89]. 
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ing the field in the late 60’s and early 70’s. Holland was the first to analyze applied 
evolutionary algorithms, and proposed what has become known as the fundamen- 
tal theorem of genetic algorithms—the schema theorem—as the principle underlying 
their success [43]. In the presence of fitness proportionate selection (to be discussed 


in Section 4.3.3), the schema theorem may be written as follows: 


m(h,t+1) > ma tE ] — po) — pmO(h) (4.1) 

where 

m is the number of times a schema (h) appears in the population, 

t is the generation index (time, essentially), 

f(h) is the average fitness of the individuals which contain schema h, 

f is the average fitness of the entire population, 

De is the crossover probability,’ 

ó(h) is the defining length of the schema in number of bits, 

l is the string length in number of bits, 

Pm is the mutation probability, and 


O(h) is the order of the schema. 


The inequality places a lower limit on the number of copies of a schema which can 
be expected to propagate to the next generation. The manner in which this limit 
is exploited by a genetic algorithm is not immediately obvious; it is first necessary 
to understand what is meant by a schema and to become acquainted with some 
important schema properties. 

A schema (plural schemata) is a similarity template common to a set of bit strings 


(i.e., encoded input or design parameters) which are identical at certain locations.? 


The schema alphabet consists of all the characters of the encoding alphabet plus a 
“don’t care” symbol. For example, the four binary strings {1000, 1010, 1100, 1110} 
can all be represented by the binary schema [1##0]. This also happens to be the 
most specific schema which represents this entire set; the remaining representative 
schemata—all of them less specific—are [l###], [###0] and [####]. These last 


three schemata represent other strings as well as the set of four above. 


T Crossover is the usual method of mixing the parameters of selected parent strings in a GA. Crossover probability 
is the likelihood that mixing will occur during a reproduction. Crossover and crossover probability are discussed 
further in Sections 4.3.4 and 4.3.5 respectively. 

8 The most common definition of a generation in a GA is a number of evaluations equal to the size of the population. 
Thus if P is the population size and the population immediately after initialization is called generation zero, then 
generation one is the population after P evaluations, generation two is the population after 2P evaluations, and so 
on. For non-generational GAs, such as the ones developed here, t may be assumed to increment each time fitnesses 
are recalculated for the population. Section 4.4 deals with generational issues in more detail. 

Some authors, such as Schaffer [78], refer to schemata as “hyperplanes”; Goldberg [33] calls them “building 
blocks.” 
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A schema can span any fraction of total string length and, as just seen, can have 
varying degrees of specificity. The schema properties which quantify these concepts 
are defining length and order. The defining length of a schema is the inclusive number 
of bits between the first and last non-“#” characters. Thus the defining length of 
schema [1##0] is four, whereas the defining length of the schema [## 101###] is 
three. The order of a schema is simply the number of fixed positions, or non-“#” 
characters, that it contains. The order of a binary schema may be determined by 
simply counting the number of 1’s and 0’s. 

Genetic algorithms proceed by selecting and copying the strings representing mem- 
bers of its population (with bias toward better individuals), swapping substrings 
among them, and occasionally mutating a bit to preserve diversity—this is the repro- 
ductive aspect. In doing so, they essentially process schemata. GA theory maintains 
that the availability of a greater number of schemata corresponds to greater efficiency 
of the algorithm’s search. The following discussion will show why this is thought to be 
true, and how a population’s schemata content can vary depending on the encoding 
alphabet. 

Simple analysis reveals that if a binary schema contains k “don’t care” symbols, 
it is represented by 2* specific strings; a higher-cardinality g-ary (q > 2) schema is 
represented by q* specific strings. Conversely, 3° schemata represent any @ arbitrary 
binary bits, although any particular sequence of £ binary bits will be represented by 
only 2° different binary schemata. For a g-ary system of the same capacity, the string 
will be £, < £ bits long, allowing the formation of (q + 1)% schemata; any particular 
sequence of l, q-ary bits will be represented by 2% schemata. For example, 3° = 27 
schemata can be formed from three arbitrary binary bits as shown in Table 4.1. Of 

000 007% 070 #00 O## #04 ##0 #FF 
001 071 #01 ##1 

010 01% #10 #1 

011 #11 

100 10# 1#0 1# # 

101 1#1 


110 11# 
111 


Table 4.1: The possible schemata of a three-bit binary string 


these schemata, only 2% = 8 are represented by any particular string. The string 
000, for example, is representative only of the schemata in the top row of the table. 
Contrast this to a coding alphabet of higher cardinality; octal, for instance. Only one 


octal bit is necessary to obtain the same capacity (eight possible values, that is) as 
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the three-bit binary string. Thus q = 8 and £, = 1 for the equivalent octal system; 
there are only nine possible schemata, these being any of the single digits 0-7 or the 
character +. Of these nine, only 2? = 2 are represented by any particular one-bit 
octal “sequence.” 1° 

In general, a bit string encoded using a low cardinality alphabet will be representa- 
tive of a greater number of schemata than the same information coded with a higher 
cardinality alphabet. A simple example of function minimization, adapted from Gold- 


berg [34] and shown in Table 4.2, will illustrate why this is advantageous. A small 


Binary Octal Function value 


000 0 8 
011 3 22 
101 5 3 
iit T 11 


Table 4.2: Binary vs. octal structures. A function value of zero is optimum. 


sample population has been evaluated by a hypothetical cost function. The function 
(objective) values—the result of sending the design parameters to the evaluator—are 
shown in the right hand column. The actual values of the design parameters them- 
selves are irrelevant to this example, but their notional binary and octal encodings 
are shown in the first two columns. A scan of the octal encodings and their function 
values reveals no obvious correlation; it is difficult to imagine how one would proceed 
with the search based on this information. The binary strings, on the other hand, 
offer more insight. It is possible that the first string is highly fit because of the 00 in 
its last two positions or perhaps because of the rightmost 0 alone. Alternatively (or 
perhaps in addition), it may be that the first and third strings are both relatively fit 
because they have a 0 at the middle position. Whatever the actual correlation may 
be, it is apparent that the correlation possibilities are increased when the cardinality 
of the encoding alphabet is decreased. This observation holds in general. The amount 
of information contained in a lower-cardinality encoding is in effect greater than that 
contained in higher-cardinality encoding, even though the non-encoded parameters 


themselves are the same in each case. In Goldberg’s words, 


. many hypotheses can be formulated regarding the association between 
substring values and high fitness [when low cardinality alphabets are used], 


10 There are obvious problems in attempting to make such a comparison between arbitrary bases; base 4 and base 5, 
for example. Both require at least two bits to encode eight values but are under-utilized to different extents (base 4 
has the capacity to encode 16 values in two bits; base 5 will allow 25). Suffice it to say that higher cardinality in 
general means fewer representative schemata, and that an under-utilized alphabet presents implementation as well as 
analysis problems. 
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and it is this information that is recombined to speculate on possibly better 


structures during the normal course of genetic search. [32] 


By evaluating individuals in the population, a genetic algorithm is in a sense 
sampling all of the schemata present in their bit strings. The greater the number of 
schemata in those strings, the more information the algorithm acquires per evaluation. 
The significance of the schema theorem now becomes more clear. Schemata which 
tend to be present in above-average individuals will promulgate and multiply while 
those which tend to be present in below-average individuals will diminish. 

The rate of propagation is also affected by the middle term in the brackets of Equa- 
tion (4.1), the probability of schema disruption by the crossover operator (crossover is 
the mechanism by which substrings from selected parents are mixed during reproduc- 
tion). This term favors schemata of short defining length ô, but makes no distinction 
among encoding cardinalities—the numerator and denominator of this term decrease 
in proportion, for a given schema, as cardinality decreases. The final term in the 
inequality, the mutation operator, further reduces the propagation rate of otherwise 
fit schemata but is somewhat of a necessary evil. It is intended to allow all possible 
schemata some chance of being produced and evaluated, even those that do not exist 
in the initial population and cannot be manufactured by the crossover operator. 

If the more fit individuals in the population are given preference in breeding and 
the less fit are more likely to die—a situation realizable with some sort of fitness- 
biased selection process—the schema theorem dictates an exponentially increasing 
number of trials for short, low-order, above-average schemata as the generations pro- 
ceed. Holland conservatively estimated this rate to be O(n*), where n is the number 
of function evaluations [32]. This gives the GA tremendous processing “leverage” 
compared to other search mechanisms; Holland named this leverage implicit paral- 
lelism. The availability of a large number of favorable schemata multipying at such a 


rate is obviously conducive to the algorithm’s search efficiency. As Schaffer puts it, 


. [implicit parallelism] constitutes the only known example of combina- 


torial explosion working to advantage instead of disadvantage. [78] 


4.3 GA Functions 


The functions of a typical GA are surprisingly simple. In fact, it is partly this simplic- 
ity combined with their remarkable search power that makes GAs attractive. Sim- 
plicity also makes for relatively easy programming and thus invites variation, as a 


custom GA can realistically be written for any new application that comes along. 
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Most original GAs—even if they employ canonical basic fuctions—have their own 
quirks and novelties, some due to their particular application and some based on 
the intuition of the programmer. GAs in the literature thus share basic features but 
are usually unique in some way, more often in implementation than in theory. This 
makes it difficult to define what exactly is meant by a “typical” GA. Nonetheless, the 
concept of a typical GA is needed here so that it can be compared to the Pareto GA 
developed during this research. For this purpose, a typical GA will be taken to mean 
one having all of the most common features among those in the current literature. 

A typical GA minimizes a single objective. The evaluator or objective function re- 
turns a measure of “goodness” for each individual and the algorithm seeks to increase 
the goodness of its population through recombination (swapping of substrings, with 
mutation). The less common multi-objective GA (MOGA) has a notable difference: 
the evaluator returns more than one objective value per individual; this objective 
“vector” is then scalarized to result in a single measure of goodness. Pareto GAs, the 
least common of the three types, also differ from typical GAs primarily in the way 
they measure goodness.!! In contrast to scalarized MOGAs, however, goodness in a 
Pareto GA is measured not by combined objective values but by relative dominance. 
To a Pareto GA, all non-dominated individuals are equally good, regardless of their 
absolute objective values. 

Apart from the number and assessment of objectives, basic functions of typical 
and Pareto GAs are essentially the same. These are shown in the flow diagram of 
Figure 4-1 and discussed in the following sections. In most cases, a brief description 
of the function will be given, followed by a mention of how it might be implemented 
in a typical GA. Any differences between the typical implementation and that used 
in this research will also be noted. Application-specific issues resulting from the use 


of DPLL as the evaluator will be covered in detail. 


4.3.1 Establishing Generation Zero 


Design Parameters, Ranges and Discretization 


An initial population, or generation zero, is usually established by generating ran- 
dom combinations of input (design) parameters, evaluating them and storing the 
results until the desired population size is attained. During the random generation 
of each new individual, the input parameters are allowed to take on any one of a 
pre-determined set of values with equal probability. A full set of input parameters 
defines a variant, or individual, and the evaluation of a set of parameters determines 


1! Schaffer [76, 77] was apparently the first to adapt a GA to search for Pareto optima in 1984. 
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establish initial population 
parent selection 


reproduction 





Figure 4-1: Flow diagram for a typical genetic algorithm. 


its objective value (or values, in multi-objective optimization) and in some cases its 
feasibility. 

For example, if one sets out to design a crate which maximizes internal volume 
while minimizing surface area, the likely design parameters are the crate’s length, 


width and height. The set 
eprom — 4 | 
defines one particular variant which, when processed by the evaluator functions 


volume = lx wxh 


surface area = 2[(1 x w) + (Lx h)+ (wx h)}, 


yields a volume of 24 and a surface area of 52. If the design requirements specify 
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a volume of, say, at least 18 and a surface area of no more than 50, evaluation has 
revealed the objective infeasibility of this variant as well as its objective values. In 
some cases feasibility may be a function of design parameters as well as objective val- 
ues; for example, the crate may be required to fit through a door of given dimensions. 
This would place upper limits on the crate’s length, width and height, effectively 
establishing their feasible ranges since zero is the obvious minimum for all three. In 
cases where no obvious input limits exist, they must be determined by heuristics or 
simply guessed and adjusted later if necessary. Once established, parameter ranges 
must be discretized such that only particular values within each range are allowed. 
Allowable values of design parameters must be discrete and pre-determined be- 
cause of the mechanics of the GA reproductive process. During reproduction, input 
parameters are converted to ordinals and then encoded into bit strings. The result- 
ing strings are commonly referred to as chromosomes, and all the bits making up a 
single parameter value are collectively known as a gene.!* Table 4.3 shows how the 
allowable values of a notional real design parameter (the length of the crate in the 


previous example, say) might be represented as ordinals and as binary genes. Width 


Parameter Ordinal Binary 


2.015 0 00 
3.030 l 01 
4.045 2 10 
5.060 3 11 


Table 4.3: Notional real parameters, as ordinals and binary genes 


and height would be encoded similarly (but not necessarily with the same resolution) 
and the resulting genes strung together to form a chromosome. If each design param- 
eter were not restricted to discrete values within a given range (for example, if the 
length of the crate were allowed to take on any real value between 2.015 and 5.060), 
the number of bits required for the gene could not be determined beforehand.'* Genes 
(and therefore chromosomes) of unspecified length present several difficulties, partic- 
ularly in mating, and make the programming of the GA much more cumbersome. 
Although such problems are not insurmountable (“messy” GAs have been developed 
which successfully operate on chromosomes of variable length [37]), the documented 
benefits of variable string lengths are not conclusive. This research will deal only with 
" 12Much of GA terminology is understandably borrowed from the biological sciences. Two less frequently used but 
related allegorical terms are allele, meaning the actual value of a bit, and locus, meaning a bit’s position on the 
chromosome. 

131¢ should be noted that discretization of parameters does not place any limitations on resolution beyond the 


precision and capacity of the hardware used. Greater resolution simply requires longer chromosomes and generally 
slows the convergence of the algorithm, as is typical of any numerical method. 
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chromosomes of equal length. 

While the feasible ranges of input parameters may be known or estimated in most 
cases, little guidance is generally available concerning their proper resolution. Obvi- 
ously the choice of base is a factor; in binary coding one would be well advised to 
choose some power of two for each parameter’s resolution. Beyond this, the choice 
must be based on the trade-off between resolution and computation time. Discretiz- 
ing each parameter to maximum precision will likely result in very long chromosomes 
and search times. GAs are very efficient search engines, but they cannot be expected 
to sample enormous design spaces effectively in limited time. When determining pa- 
rameter discretizations, then, it is advisable to keep the time vs. precision trade-off in 
mind. Lengthening a binary gene by one bit, for example, doubles the size of the de- 
sign space. Discretization decisions should generally favor low to moderate resolution, 


at least until preliminary results have been reviewed. 


DPLL Parameters 


Table 4.4 shows the design parameters used in this research with DPLL as the eval- 


uator; they are described in detail below. Parameter ranges here are chosen via a 







Input Parameters 


COo o np Parami O O 
[configuration | 14 f 4 f meer 
“tip radius | 0.02:0.05 | 8 |% body length 
duct/rotor intersect | 02:08 | 8 | %educt chord 
[duct chord length | 0.001 :0.080 | 8 [Y body length | 
duct load Sid 004.0 | 8 [| 6 — 


Table 4.4: Input parameters and ranges used for DPLL optimization. Blade circulation values are 
positive for rotors, negative for stators. 










combination of heuristics and trial runs so as to allow all reasonable possibilities. Most 
parameters are given a resolution of eight; this is assumed to be a reasonable initial 
compromise between precision and computation time. The greater resolution of the 
fullness factor is due to its strong correlation to one of the objective values (volume); 
this will be discussed further in Chapter 5. For the purposes of this optimization, 
other inputs required by the stand-alone version of DPLL (forward speed and number 
of propeller blades, for example) are either fixed for all variants or are calculated as 
functions of the parameters shown. The total number of variations possible (i.e., the 


size of the design space) using this resolution is 1.77 x 10%, 
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1. Configuration determines the number, type and relative locations of propeller 
stages. Configuration 1 is rotor-only, 2 indicates a stator upstream of the rotor, 
3 astator downstream of the rotor, and 4 indicates two stators, one upstream and 
the other downstream of the rotor. These four different configurations are treated 
by the reproductive mechanisms of the GA as separate, co-located species. They 
are not allowed to interbreed but do compete for the same “resources” (fitness 
values) and are therefore subject to dislocation and extinction. The concept of 


competitive species is discussed further in Section 4.4. 


2. Fullness factor determines the shape of the stern. Each variant begins with 
a common set of B-spline vertices defining the body profile. This template 
produces a gently sloping profile from mid-body to stern and is not modified if 
the fullness factor is 0.5 (the minimum). Any fullness factor greater than 0.5 
scales up the z coordinates of the vertices defining the stern, but leaves the r 
coordinates unchanged.!° The result is a greater extent of parallel mid-body and 
a more abrupt transition at the stern. The relationship between vertex scaling 
and fullness factor is set such that a fullness factor of 1.0 (the maximum) results 
in a transition abrupt enough to cause separated flow for nearly all combinations 
of the remaining parameters; this is one instance where trial and error is required 
in setting parameter limits. Figure 4-2 shows the range of hull profiles generated 
by varying the fullness factor between 0.5 and 1.0. The actual resolution is 
about twice that shown, as intermediate profiles have been omitted from the 


plot to avoid clutter. 


3. Tip radius is the r coordinate, as a fraction of body length, of the governing 
propeller tip. It also specifies the radial offset of the duct, as the duct camber line 
is required to intersect the governing propeller tip (DPLL uses a preliminary rigid 
translation of the input camber line to ensure this intersection). In the modified 
version of DPLL used in this optimization, the governing stage is automatically 
specified as the one which is farthest forward on the body. Tip radii of any 
additional stages are calculations based on the solved shape of the duct camber 


line during each iteration. 


4. Hub circulation specifies the non-dimensional local circulation on the lifting seg- 
ment nearest the hub. As discussed in Section 2.2.1, circulation is essentially the 


14 There is no particular reason for the lower limit being 0.5 other than to avoid a lower limit of zero, which might 
be construed as zero volume. In any regard, the numerical value itself is meaningless; it simply represents one of eight 
possible sets of body profile vertices to DPLL and an ordinal number, for encoding, to the GA. 

!5 The vertex which closes the stern at z = 1.0 and the vertices defining the sting are unaffected (the sting is 
described in Section 2.3.7). 
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Figure 4-2: Range of allowed hull profiles for optimization. 
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velocity difference across a foil and is thus a measure of blade “load.” Indepen- 
dent hub circulation values are specified for all three propeller stages, regardless 
of how many are actually present (that is, regardless of configuration number), 
in order to maintain uniform chromosome length. Extraneous values are simply 
disregarded during evaluation. Circulation values for stators are made negative 


so as to induce tangential velocities opposing those of the rotor. 


. Tip circulation is the non-dimensional local circulation on each stage’s outermost 
lifting segment, specified in the same manner as the hub circulation. Values 
of circulation at intermediate lifting segments are determined by interpolating 
between the hub and tip values, resulting in a linear distribution across the blade 


span. 


. Duct/rotor intersect specifies the intersection of the duct camber line and the 
rotor tip as a percentage of duct chord length. It effectively positions the duct 


axially, as the rotor location is fixed and constant for all variants.!® 


. Duct chord is the axial distance between the duct leading and trailing edge, as 
a fraction of body length. If the configuration value is 1 (no stators present), 
the minimum duct chord allowed is 0.001L. This negligible value is used to 
approximate a non-ducted system, as the current version of DPLL requires the 
presence of a duct. For configurations which include stators, the no-duct ap- 
proximation is not allowed. This is because all propeller stages must be enclosed 
by the duct, and negligible duct chord would result in an unrealistically small 


separation between successive stages. 


. Duct load is non-dimensional total circulation on the duct, defined as the sum 
of the camber line vortex ring strengths. This definition of duct load is closely 
correlated to lift; a negative load with a positive angle of attack will produce 


negative thrust (i.e., opposing the rotor). 


Sample Chromosome 


Table 4.5 shows a notional binary chromosome, such as might be generated by the 


initialization or reproduction routines and sent, in its real form, to DPLL for evalu- 


ation. The configuration gene (cfg) of this variant decodes to ordinal 2, meaning 


that it is the third of the four possible configurations and that its propulsor consists 


16 Rotor location is an example of a design parameter which is fixed simply to keep combinatorics to a manageable 
level. It would no doubt be an interesting parameter to vary in further studies. 
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Hsipo 
CI O E aT 





A 3 
real [3.0 [ 0.7100 |0.0286/0.0135/0.0207|0.0314|0.0243/0.0171 


Table 4.5: Example chromosome with decoded values. 


of a rotor followed by a stator.!” 

The fullness factor (fulfac) decodes to ordinal 7, out of 16 possible values. The 
real number, an arbitrary scaling of the ordinal to the range 0.5-1.0, is 0.71 by 
interpolation. 

The tip radius decodes to ordinal 2, out of 8 possible values between 0.022 and 
0.05Z inclusive (LZ is body length). Thus the real tip radius by interpolation is 
0.0286L, or 2.86% of body length. Since this is a rotor-stator configuration, the 
rotor is the most forward and therefore the governing propeller stage to which this 
tip radius applies. 

The first and second hub and tip circulation values are assigned to the rotor and 
stator respectively; the third hub and tip circulations in this chromosome are ignored. 
The rotor tip circulation gene (Hip) will be used here as an example. This gene 
decodes to ordinal 4, out of 8 possible values between 0.010 and 0.035 inclusive. The 
real non-dimensional circulation A is thus 0.0243 by interpolation. Using the decoded 
value of the tip radius gene from above and Equation (2.12), the local tangential 
velocity at the tip is found to be 0.106V, or about 10% of ship speed (as for the 
remaining unknowns in Equation (2.12), all variants considered in this optimization 
have a maximum hull radius of 0.052 and five propeller blades). Note that this is 
induced velocity only and does not include the rotor’s angular velocity. The remaining 
hub and tip circulation genes decode similarly, except that the radial coordinates of 
these points vary slightly as the iterations progress and the solved shape of the duct 
camber line changes. 

The duct/rotor intersect gene (Syotor) specifies the axial position of the duct camber 
line so that it intersects the rotor tip at 45.7% of camber line arclength. The duct 
chord gene (Cguct) indicates that the chord length is 0.0461Z, or 4.6% of body length. 
Since the rotor location is fixed, these two parameters along with the tip radius fully 
describe the initial position of the duct. Stators, which must be individually located 
with respect to the duct chord when DPLL is run as a stand-alone program, are 
automatically positioned in this optimization. For the rotor-stator configuration of 
MrT helo possible cenfeurations dictates two-bit binary gene, which has ordinal values of zero through three. 


Configuration numbers are the ordinal values plus one; this is simply an aesthetic avoidance of zero as a configuration 
number. In retrospect, the translation is probably more confusing than helpful. 
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this example, the stator tip is placed midway between the rotor and the duct trailing 
edge and the lifting line is grown inward, normal to the local velocity, to intersect the 
hull profile. 

The duct load gene Gguce decodes to a non-dimensional circulation of 0.06, which 
is essentially the velocity delta between the duct’s upper and lower surface, non- 
dimensionalized by the radii of the duct control points and ship speed and summed 
over all duct control points. The positive value indicates counter-clockwise circulation 
about a centerline (9 = 0) section through the duct, so by Equation (1.1) it should 


produce lift with a component in the forward (—z) direction. 


Population Sizing 


Population size is among the parameters having the greatest impact on GA perfor- 
mance, second perhaps only to selection method. There are no firmly established 
guidelines for determining the correct population size; in fact, sizing is somewhat of 
a dilemma. If the population consists of only a few individuals, the selection process 
is hindered by a lack of choices. The number of schemata available to the GA is 
very limited and improvement will almost certainly be intolerably slow. If by some 
fortunate accident an exceptionally fit individual is produced, it stands a fair chance 
of being quickly eliminated due to the quirks of stochasticity. At the other extreme, 
if the population size is huge (approaching the size of the design space, say), then 
the GA is irrelevant because establishing generation zero requires a near-exhaustive 
search and the optima will become evident as by-products of this process. 
Population size also affects convergence time. Goldberg [35] has shown that typical 
GAs converge in O(log P) generations, where P is the population size and a generation 
involves P cost function evaluations. This would indicate that smaller populations 
require fewer evaluations. Convergence, however, is not necessarily equivalent to 
locating the global optima. Again considering the extreme cases, a two-member 
population will probably converge (stop improving, that is) quickly since there are 
very few schemata to process. Convergence will most likely occur at a point near the 
original two members. In contrast, a huge population will be far from converging after 
the same number of evaluations, but will probably have better solutions available at 
that time due to the extensive initial “random search” performed while establishing 
the population. These general observations are supported by the work of De Jong, 


who performed some of the first experiments on GA convergence properties in 1975: 


Increasing the population size was shown to reduce the stochastic effects 


[of random sampling on a finite population] and improve long-term perfor- 
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mance at the expense of slower initial response. [19] 


The “correct” population size in terms of efficiency and efficacy obviously lies 
somewhere between the two extremes. De Jong suggested a population size of 50- 
100, based on his suite of test functions. Grefenstette [40] took the novel approach 
of using a GA to determine optimal GA performance parameters and suggested a 
population size of only 30. In one of the first analytical investigations, Goldberg [31] 
proposed an optimal population size based on the length £, in number of bits, of the 


chromosome: 
O 2 same (4.2) 


Schaffer et. al [79] later proposed an empirically optimal P of 20-30, and Goldberg 
in a more recent analytical work suggested that with a few simplifying assumptions, 
the correct population size is probably of O(£) [36]. 

A population size of 200 is used in this research. The chromosome length is 
36 binary bits, which gives P = 311 using Equation (4.2). All of the remaining 
studies mentioned above, including Goldberg’s more recent work, would indicate a 
much smaller population size; however, the Pareto GA developed during this research 
operates on four distinct species simultaneously (see Section 4.4). As this is an original 
concept, there is no way of predicting its effect on optimal population size. The value 
of 200 is assumed to be conservatively high; paying an initial response time penalty 


is considered preferable here to premature or false convergence. 


4.3.2 Evaluation 


Evaluation is the process of determining the objective value(s) of a member of the 
population, represented as a set of design or input parameters. Individuals are sent 
to the evaluator as a set of inputs in evaluator language, meaning that they must be 
converted from their encoded chromosome states back into actual parameter values. 
The evaluator is completely independent of the other GA functions; it has no memory 
or biasing features of its own. It is often referred to as the objective or cost function, 
although it need not be an explicit function in the mathematical sense and often is not. 
In fact, it can take almost any form, such as complex stand-alone analysis programs 
(such as DPLL) or even experimental testing (for an example of this, see Barrett’s 
application to flexible underwater vehicles [2]). Objective values returned by the 
evaluator provide exclusive guidance for the GA’s search process. This is in contrast 


to the majority of optimization algorithms, which utilize gradient information. 
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In this research, DPLL is used as the evaluator. The input parameters have been 
previously discussed; the calculated values of interest—the objective values—are us- 
able hull volume, power coefficient and, eventually, cavitation index. Hull volume 
is determined by a trapezoidal integration of the body profile and adjusted by es- 
timates of internal reduction gear and shafting volumes as functions of developed 
torque. Appendix B details the method used to estimate these internal volumes. 
Power coefficient is defined as: 

Qu 


Ch = 
—— JpV3Ap 


(4.3) 


where Q is shaft torque, w is shaft angular velocity, V, is forward speed of the ve- 
hicle and Ag is the cross-sectional area of the vehicle at the point of maximum 
radius. Shaft torque and angular velocity required for force-balancing are calculated 
by DPLL; forward speed and the body profile (and thus the maximum cross-section) 
are inputs. Since forward speed and maximum radius are identical for all variants 
in this optimization, the power coefficient serves as a measure of relative propulsive 
efficiency. Cavitation index is a function of lift at the rotor tips, relative fluid velocity, 
tip speed, and geometry. Appendix A describes this calculation in detail. A lower 
value of cavitation index as defined there corresponds to reduced cavitation and noise. 

Thus the general objectives of this Pareto optimization—not to be confused with 
the objective values of any particular variant as per DPLL—are an increase in vol- 
ume, a reduction in power and a reduction in cavitation at the given forward speed. 
Superlatives such as maximum or minimum are inappropriate, as the achievable lim- 
its are initially unknown and the GA operates only on comparisons. Even if the 
limits were known, the GA could not guarantee discovery of the corresponding design 
parameters, being non-deterministic. This is why GAs in general have no inherent 
termination point. Among the typical convergence criteria used are some number of 
generations without improvement (this corresponds to frontier stability in Pareto GAs 
and is the criterion used here), the attainment of some pre-defined objective goal, or 


) 


the onset of “population takeover,” where the population becomes homogeneous or 
nearly so. This last phenomenon occurs only in GAs which allow the existence of 


duplicates (more about this later). 


4.3.3 Selection 


Selection is the process by which members of the population are chosen to mate and 
reproduce. Selection and breeding of individuals with above-average schemata, as 


evidenced by their superior objective values or relative dominance, should generally 


100 


produce above-average progeny, thus moving the population as a whole toward im- 
provement. The most obvious method of selecting parents, simply choosing the best 
phenotypes, produces poor results because the search quickly becomes localized.!® 
This often leads to premature convergence near the best members of the initial pop- 
ulation. In order to achieve the previously mentioned balance between exploitation 
and exploration, a selection method must favor the more fit individuals—this requires 


biased-random selection based on fitness. 


Objective Values and Fitness 


Objective value and fitness are generally not equivalent. An individual’s objective 
value is a function only of its design parameters, with the function being defined 
by the evaluator. Fitness, on the other hand, is a relative term which indicates an 
individual’s standing among its contemporaries. More formally, if Y is a set of input 
parameters and O = E (b) is the vector of objective values returned when wp is sent 


to the evaluator, then the fitness F of Y is given by: 
FẸ) = f (E), t) = £(6,2) (4.4) 
where t is the generation index of the algorithm. When all individuals in a population 


have been assigned fitness values, selection proceeds by making stochastic choices 


among them based on fitness. Those selected go on to produce offspring. 


Fitness Scaling, Clones and Ranking 


One of the simplest fitness functions f to implement is a scaling of raw objective val- 
ues to the limits of the current population. This is known as proportionate or linear 
fitness scaling [41]. Fitness scaling followed by some form of stochastic selection is 
a decided improvement over simply selecting the individual with the best objective 
value; however, the problems of intense local selective pressure and premature conver- 
gence are not entirely eliminated. If one member of the population happens to be far 
superior to the others, fitness scaling tends to favor that individual to the exclusion 
of others. The result may be a population full of duplicates or near-duplicates of one 
early, fortuitous (but non-optimal) individual. 

Some discussion is in order here about duplicates, or clones, and how they affect 
population takeover. Clones are individuals which are identical in terms of all ac- 


tive design parameters to some other member of the population. They usually arise 


18 Another example of borrowed biological jargon. The GA genotype is the set of design parameters (the evaluator 
input); usually it refers to an encoded chromosome but may also mean the set of non-encoded parameters. The 
phenotype is the set of corresponding objective values (the results of an evaluation). In biology, the genotype is 
the genetic code (an organism’s chromosomes, basically) and the phenotype is the manifestation of the code (the 
expression of traits, or the characteristics of the organism) [85]. 
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through reproduction but may occur during the initialization process as well. In 
single objective GAs and MOGAs, they do not present a major problem as long as 
the population size is sufficiently large and some mechanism, such as mutation, is in 
place to preserve diversity and prevent premature convergence. The proliferation of 
clones to the point where they comprise the majority of the population, making fur- 
ther improvement unlikely, is known as population takeover. In single objective GAs, 
this is expected behavior and is often used as an indicator of convergence. In Pareto 
GAs, however, clones are usually considered detrimental because they decrease the 
potential frontier resolution and population diversity [91]. 

There are three basic ways of dealing with duplicates: they can be tolerated, penal- 
ized or rejected. Since tolerating clones is counter-productive in Pareto optimization, 
the choice here must be between penalization and rejection, even though both meth- 
ods require additional processing time for clone identification.!? The only rationale 
for choosing penalization is that time already spent generating the clones is wasted 
if they are then identified and rejected.*° This argument is not persuasive, consid- 
ering that (1) even more time will be spent determining and assigning the penalty, 
(2) the problem of diminished resolution is only partially solved since some number 
of clones will generally be present in the population, and (3) reproduction, mutation 
and identification time—that is, the time required to generate and identify a clone— 
is usually insignificant compared to evaluation time. For the Pareto GA developed 
here, the evaluation of a single individual by DPLL requires at least ten minutes on a 
400 MHz DEC Alpha workstation, while all GA functions combined (including clone 
identification and ranking of the entire population) require less than two seconds. 
Penalization is therefore not justified in terms of time savings, and rejection is used 
exclusively here. 

Returning now to the discussion of proportionate fitness assignment: recall that 
proportionate assignment followed by a stochastic selection method alleviates some 
convergence pressure, but fitness remains based on relative phenotypical distance. 
This is not conducive to exploration and can lead to poor results, particularly for 
multi-modal functions. Even single objective/single optimum searches can be mislead 
by fitness scaling if a relatively fit individual arises in a flat, non-optimal region of 
objective space. Also, since fitness scaling requires a scalar objective value, it must 
be preceded by scalarization (a weighting method, for example) when used on multi- 
objective problems. It is therefore not compatible with a true Pareto GA, and is not 
~ 19§chott (80] performed convergence studies on a Pareto GA using each of these three methods, confirming that 
both penalization and rejection are significantly more efficient than tolerating clones. 


20 Note that this refers to reproduction, mutation, or initialization time, not evaluation time. There is no conceivable 
reason for evaluating an identified clone—its objective values are already known. 
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suitable for this research. 

One way of smoothing out the intense pressures of fitness scaling is to employ 
a ranking scheme. This involves assigning fitness according to relative standing in 
the population, regardless of the phenotypical distances involved (e.g., best = 1, 
second best = E, W. worst = pr If two or more individuals have equal standing, 
they are assigned an averaged rank. For example, if there were n “best” members, 


they would all be assigned a rank of 


P+(P—1)+(P-2)+-:-+(P—n41) _ | =o 
+ 51-28 (4.5) 


Stochastic selection then proceeds as it would with proportional fitness scaling, but is 
based on rank rather than relative objective value. Ranking methods are well suited 
for Pareto optimization, where relative dominance is the only meaningful criterion for 
comparing individuals. The “best” members of the population in the Pareto sense 
are those that are non-dominated, the “second best” are those dominated only by the 
“best,” and so on. 

Ranking also allows control over the ambient selection pressure [91]. For example, if 
the simple ranking system mentioned above provides insufficient pressure—as might 
be indicated by excessive convergence times—the assignments can be made on an 
exponential rather than a linear ranking scale. The general implementation would 
be to assign fitnesses of 1, s’, s’t!, ..., from best to worst, where s < 1.0 and j is 
some positive integer. Both s and j may be varied to obtain the desired level of 
pressure [27, 41]. Ranking methods employed in this research use such a scale, with 


P-1 : 
s = ao and J] = 3. 
Roulette Wheel and Tournament Selection 


Stochastic selection mechanisms following either a proportionate fitness assignment 
or a ranking scheme appear in several variations in the literature. The most simple 
of these, used in this research, is known as “roulette wheel” selection. A vector of 
cumulative fitness values for the population, in arbitrary order, is constructed as 


follows: 


P 
F= Cr i): Sai a nt (4.6) 
ere hi 
A random number r of uniform probability distribution between zero and one exclu- 
sive is then generated and the member whose index in F is defined as follows is the 


chosen parent: 


ME JF eee with y =0 (4.7) 
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A drawback of the roulette wheel method is that it is inherently noisy; 1.e., the 
variance of the probability density function for selecting a string of given fitness 
is relatively large. This noise, along with that generated by the reproduction and 
mutation operators, tends to obscure the signal difference between good and bad 
schemata, thus hampering the GA’s ability to select well. Goldberg [36] has analyzed 
the effects of noise and compared some of the more common selection methods using 
a suite of test functions of varying complexity. His results, in addition to confirming 
the advantages of ranking over fitness scaling, also show that some alternatives to 
roulette selection may be preferable. 

One such alternative is tournament selection. A tournament is held among some 
pre-defined number of individuals—chosen randomly from the population—and the 
tournament winner is selected as a parent. No roulette-type routine is necessary, as 
the random selection of contestants provides the required stochasticity. Tournament 
victory can be based on scaled fitness values (which, for tournament methods, is 
equivalent to simply using raw objective values) or on relative rank. Dominance 
ranking is the only viable alternative when applying tournament selection to Pareto 


GAs. 
Niching 


Rank-based tournaments do not necessarily result in a unique winner and therefore 
must include some means of resolving a tie. Selecting randomly from the tourna- 
ment winners is reasonable, but there are advantages to be had by applying greater 
discretion. These advantages are a result of the fact that multiple variants in close 
phenotypical proximity (those which are crowded together in objective space) con- 
tribute little additional information to the GA. In fact, they can be detrimental 
to performance as they take up space in the population which could be occupied by 
more diverse variants. In this sense, such “near neighbors” are quite similar to clones, 
although phenotypical proximity does not necessarily correspond to genotypical prox- 
imity. Pareto GAs are particularly affected by phenotype crowding, as their goal is a 
uniform distribution of the population across the non-dominated frontier. 

Resolving tournament ties and discouraging phenotype crowding are compatible 
objectives. Ties can be broken by selecting the individual from among the winners 
whose region in phenotype space is the least crowded. The degree of crowding around 
an individual is determined by conducting a “niche count.” The niche count is simply 
the number of other members which lie within some hyper-volume, or niche, centered 


on the individual in objective space. The obvious parameters necessary to define such 
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a niche are its size and shape.“ Although there has been some investigation into the 
efficacy of various niche definitions, including asymmetric shapes [45], there remains 
little evidence as to the benefit obtained by such complications. 

The simplest approach to niching—the box niche used here—utilizes pre-defined 
phenotypical distances to test for proximity. For example, individual X is said to 
be located within niching distance of individual Y if the difference in their values of 
objective Á is less than some Aniche and the difference in their values of objective B is 
less than some Bhniche- These tolerances may be static or dynamic according to latest 
population statistics. In either case, they are usually based on a calculated optimal 
separation distance, determined by dividing the assumed frontier surface area by the 
size of the population [27]. Note that if X is in the niche of Y then the converse is 
true, and that as the number of objectives increases, the time-averaged niche density 
decreases for a given population size. 

Niching should be a concern in any Pareto GA, not only those employing tour- 
nament selection. Niche counts may be accounted for in ranked GAs, for example, 
by applying the count as a rank penalty prior to stochastic selection. This is the 


approach used in the ranking methods of this research. 


Constraints 


The selection process includes the handling of constraints on objective values. The 
crate example of Section 4.3.1 demonstrated that objective feasibility can only be 
determined by evaluation. What the example’s simplistic objective functions failed 
to emphasize, however, is that significant computation time may have already been 
invested before a variant is known to be infeasible. This is certainly the case in this 
research, where a single evaluation by DPLL can require up to 30 minutes. Time 
invested must be a factor in determining how to deal with infeasible individuals. 


Three of the more common methods used are [18]: 
1. immediate rejection 
2. repair 
3. penalization 


Rejection—simply discarding the offending variant—guarantees the feasibility of the 
population, but is inefficient. Since no consideration is given to the degree of violation, 
rejection ignores the potentially useful information present in an “almost feasible” 


21See Fonseca and Fleming (26] for a description of how niche shape can be affected by varying the degree of the 
Holder metric. 
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individual. It is therefore not well-suited for identifying the limits of feasibility and 
may overlook constrained optima. Another concern is that for problems having highly 
constrained solution spaces, the identification of any feasible solutions may be nearly 
as difficult as finding the optimal solutions. Since rejection extracts zero information 
from infeasible variants, it is likely to produce poor results in such cases [33]. 

Repair of infeasible variants is necessarily problem-specific and complicated. The 
type and extent of the repair must depend on which constraint has been violated and 
the magnitude of the violation; the repair algorithm, in order to be effective, must 
be a complex, dynamic process which “learns” and improves itself as the generations 
progress. Repair is not necessarily more efficient than rejection, as the repaired vari- 
ant will require re-evaluation and the repair process itself may be computationally 
intensive. 

The penalty method assumes that infeasible variants may contain useful genetic 
information and allows them to remain in the population, although with a fitness 
penalty. This method is more likely than rejection to locate constrained optima, as 
variants which are “barely” infeasible may provide schemata which lead to the optima. 
Unfortunately, the assignment of penalties is necessarily somewhat problem-specific. 
Over-penalization will produce essentially the same results as rejection, but with an 
added computational burden. Under-penalization can be even worse: if the penalty 
is not at least equal to the objective advantage gained by violating constraints, the 
population will quickly become saturated with infeasible members. It is not clear 
how a balance between these two undesirable extremes can be struck without prior 
knowledge of the solution space or some amount of trial and error. 

This research employs the penalty method and originates a technique which may 
be helpful in generalizing penalty assignments. The constraint in this case is flow 
separation; it is violated if separation occurs upstream of the rotor. The requirement 
of attached flow may be considered a “soft” constraint, as the degree of violation 
increases with the upstream distance of separation rather than being a simple fea- 
sible/infeasible dichotomy. Variants which violate this constraint are penalized by 
decreasing their relative dominance in the population. The penalty is applied as an 
integer number of dominating “phantom” members, seen only by the individual to 
whom the penalty applies. This number increases with the distance between the ro- 
tor and the upstream separation, and has a dynamic lower limit equal to the total 
number of infeasibles in the current population. This floating lower limit is the gener- 
alization developed here and mentioned above; it essentially acts as negative feedback 


on the proportion of infeasibles in the population. If the number of infeasibles begins 
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to increase—possibly because the magnitude of the feasibility penalty has initially 
been set too low—the minimum penalty for infeasibility also increases. This gives 
a greater selective advantage to the feasible individuals, which are never considered 
dominated by constraint violators regardless of the relative objective values involved. 
The result is a decrease in the number of infeasibles in the population over time, with 


a corresponding decrease of the minimum infeasibility penalty. 


Specific Methods Investigated 


Three specific selection methods are investigated here, including two variations on 
ranking and a tournament method. All employ niching and constraint penalization. 
These will ultimately be compared in terms of their ability to locate and define the 


Pareto frontier for DPLL’s volume and power calculations in Chapter 5. 


1. Goldberg ranking”? is based on dominance layers. Non-dominated individuals 
(those in the outermost layer, or layer zero, in objective space) are assigned 
the best rank, individuals which are dominated only by those in layer zero are 
assigned the next best rank, and so on. Niche counts are then performed and the 
population is again ranked, this time according to niche count. This results in two 
ranks per individual; these values are added and a penalty is applied to variants 
which violate the separation constraint. The population is sorted according to 
this combined and penalized rank. Fitness values are then assigned to the sorted 
population according to the exponential scale discussed in Section 4.3.3. Parents 


are selected by consecutive applications of a roulette wheel routine.* 


2. Fonseca and Fleming ranking assigns each individual a rank based on the num- 
ber of other individuals by which it is dominated. As with the Goldberg method 
above, these ranks are modified by niche counts and constraint violation; the 
population is then sorted and assigned fitness according to an exponential scale. 


Parents are selected by the same roulette wheel routine. 


3. Tournament selection randomly chooses 10% of the population and selects the 
dominant member. A tournament contestant which violates constraints will in- 
cur a relative dominance penalty. This penalty will reduce a constraint violator’s 
tournament standing by at least the number of infeasible contestants. The net 
effect is that a constraint violator tied for first place prior to penalty assessment 


will instead be eliminated from consideration. A constraint violator which oth- 


22 Ranking methods are named here according to the authors who proposed them. This terminology is not universal. 
23 Selection of the second parent is somewhat complicated by the presence of other species, as discussed in Section 4.4. 
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erwise would have won the tournament outright will finish in a tie for first place 


at best. Post-penalty ties are resolved by niche count. 


It should be noted that these three methods are not arbitrarily chosen from among 
several possibilities. There are apparently no other general means of applying purely 
Pareto criteria in selection—at least there are none documented [28]. The main Pareto 
GA framework developed here incorporates these selection methods as subroutines. 
Switching among them is done prior to compiling. The objective here is to main- 
tain uniformity in all non-selection functions, allowing valid comparison of relative 


performance among the three methods. 


4.3.4 Reproduction 


Reproduction is the process by which selected members of the population breed to 
produce new variants, which are not surprisingly referred to as offspring or children. 
During the selection function just discussed, highly fit individuals are chosen to be 
parents on the basis of their objective values. Since the evaluator is essentially a 
function relating input parameters to objective values, and objective values are a 
GA’s only criteria for selection, the chromosomes of selected parents will on average 
contain relatively high quality schemata. It is the job of the reproduction process to 


mix and match these schemata so as to produce more highly fit children. 


Coding 


The code to be used for converting design parameters to bit strings for mating is one 
of several controversial aspects of GA research. Holland’s work was done with binary 
coding, and this remains the most common method in the literature. However, higher 
cardinality codes have been used with success and there has recently been considerable 
interest in “real-coded” chromosomes, where the real design parameters themselves 
are used with no encoding at all (51, 94]. Any cardinality greater than two—including 
real code—complicates the implementation of the mutation operator somewhat. The 
complication is not prohibitive, however, and is offset by shorter chromosomes and a 
more rational parameter-to-gene relation. 

The fact that such methods have been known to perform well has prompted re- 
searchers to re-examine Holland’s schema theorem, which seems to favor low car- 
dinality alphabets as discussed in Section 4.2. Goldberg [34] attempted to explain 
the paradox by developing a theory of real-coded GA operation, which he called the 
theory of virtual alphabets. He concluded that real-coded GAs are successful due 


to implicit internal reduction of high-cardinality alphabets to low-cardinality virtual 
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alphabets by the selection operator. This explained the successes of real coding while 
preserving the validity of the schema theorem. In the same paper, he also showed that 
real-coded GAs are subject to an objective function-dependent phenomenon called 
blocking, which in certain cases can prevent them from locating optima that are ac- 
cessible to lower cardinality codes. In keeping with the opinions of the majority of GA 
researchers and Goldberg’s results, binary coding is used exclusively in this research. 

Even among binary codings, there are alternatives. Some authors advocate the 
use of Gray code as opposed to the more traditional base 2 [3, 11]. In the Gray code, 
any two adjacent base 10 integers differ at only one bit location. Adjacent base 10 
integers represented in base 2, on the other hand, may differ at all bit locations. 
For example, 32 = 011 and 4, = 100. Gray code is therefore thought to be more 
compatible with the mutation operator, as the “flipping” of any one Gray-coded bit 
is less likely to result in a decoded parameter far removed from the original. Slight 
performance advantages due to the use of Gray code have been reported in some 
cases [79, 44]. Regardless, the coding used in this research is exclusively base 2 in 


keeping with convention. The investigation of Gray coding is left to future research. 


Crossover 


Crossover is by far the most common mechanism for mixing parental characteristics. 
It has several variations; the simplest involves the production of two children from 
two parents using what is called “single point” crossover. A crossing, or cut, location 
is randomly selected on the chromosomes of the parents. The first child is then 
defined as the portion of the first parent’s chromosome prior to the cut concatenated 
with the portion of the second parent’s chromosome after the cut. The second child 
is defined by reversing the order of the parents’ contributions. Variations include 
increasing the number of cut points (resulting in multiple, or n-point crossover [33]) 
and “shuffle-cut-unshufHle” schemes, which are intended to allow any two genes the 
same probability of being passed together, regardless of the chromosomal distance 
separating them [70]. The latter process eliminates bias due to gene loci, which in 
most cases are arbitrary specifications by the programmer. 

The advantages, if any, of multiple point crossover have been shown to be small [79]. 
Holland’s schema theorem certainly indicates that schemata disruption should be 
minimized, and as Goldberg [36] puts it, “ ... no convincing evidence of high-order 
mixing success—em pirical or otherwise—has yet been offered.” Shuffle-cut-unshuffle, 
as with many similar proposals in the literature, is intuitively attractive but has 


little supporting data. In the interest of limiting the departures from standards in 
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this research to the areas of primary investigation (selection, competitive species and 
steady-state processing), standard single point crossover is used exclusively here. The 
only non-standard aspect, made necessary by the species concept, is to disallow a cut 
point location within or immediately following the configuration gene. Since parents 
must be of the same configuration, such a crossover point would result in useless 


clones. 


Mutation 


Mutation is the simplest of all GA functions, particularly when binary coding is used. 
When children emerge from the crossover process, but prior to decoding their bit 
strings back into real design parameters, each bit in their chromosome is subjected to 
inversion with some (normally very small) probability. This is to maintain some mea- 
sure of diversity in the population and prevent the permanent loss of any particular 
schema. In this research, the probability of mutation is set to 0.005; suggested rates 
in the literature vary from 0.001 to 0.01 [79]. Mutation is prevented from occurring 
within the configuration gene, which determines an individual’s species. Mutating 
this gene would be the equivalent of biological saltation and has no counterpart in 


natural processes. In the Darwinian sense, the species here are immutable.?* 


4.3.5 Replacement 


Typical GAs maintain a constant number of members in the population. Excep- 
tions do exist in the literature, but variable population size, like variable chromosome 
length, makes analysis and implementation much more difficult without obvious off- 
setting benefits [80]. If population size is to remain constant, offspring must replace 
existing members. In a typical GA, the replacement function is autonomous because 
the entire population is simultaneously replaced by an equal number of offspring; each 
new population becomes the next generation. The predominance of generational GAs 
in the literature is such that the term itself is used only in the few papers that propose 
alternatives [92, 80]. 

During the course of one generation, a simple generational GA of population size P 
calls the selection and reproduction routines Z times, resulting in E pairs of parents 
which produce two children each. The resulting P offspring are not immediately 
evaluated and are not considered to be members of the population until all selections 
and matings are complete, after which all new offspring are evaluated and they then 
24 Charles Dagwin.fonmulatedshisisheotwet natunasiaetangtith a a a a dae 


little to say about the evolution of full stern submarines. Nevertheless, given the nature of this thesis, it seemed a 
shame not to cite him somewhere [17]. 
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become the population. A common variation on the generational process is some form 
of “elitist” strategy, which ensures the survival of highly fit members of the former 
population at the expense of an equal number of poorly fit children. Otherwise, no 
continuity necessarily exists from one generation to the next. 

In contrast, a non-generational or incremental GA calls the selection and repro- 
duction routines p times per increment, where 1 < p < E The children are then 
evaluated and inserted into the population; the subsequent selection is performed 
based on the modified population. Obviously, a generational GA may be viewed as 
simply the special case where p = E, but there is an important distinction: For any 
p< £, the members to be replaced by the offspring (or equivalently, the members 
which are allowed to survive) must be selected. Thus incremental GAs require an 
additional function; perhaps this explains their relative scarcity. Given that the pro- 
cessing time for GA functions is usually quite short compared to that of evaluation, 
however, this requirement should not exclude incremental GAs from consideration. 
Parent selection routines, which must be present anyway, can be easily adapted to 
select unfit members for replacement, possibly even concurrently. Also, for purists, 
an incremental GA is more in keeping with natural population dynamics, as natural 
populations are not known for mass simultaneous births and deaths.* 

Intuitively, an incremental process would seem to provide higher quality feedback 
to the algorithm, as each new selection is made using the latest data available. This 
observation has been made by Whitley [91], and an incremental process is used in 
his GENITOR algorithm. Schott also applied incremental strategy in Pareto opti- 
mization of fault tolerant systems [80]. There are, of course, critics of this approach. 
Hancock [41] claims that incremental reproduction inevitably involves the same kind 
of sampling errors as roulette wheel selection, and a study of selection methods by 
Goldberg and Deb [35] concludes that the GENITOR algorithm produces very high 
selective pressure. According to Hancock, however, this pressure is primarily due to 
perpetual replacement of the worst individual in the population by new children, a 
feature which is obviously not inherent in incremental strategy. 

The GA implementations written for this research investigate the case where p = 1; 
that is, where children are immediately evaluated and inserted into the population. 
This variation on incremental strategy is known as steady-state processing. In addi- 
tion to being intuitively attractive, this method is chosen to avoid instability when 
tournament selection is combined with niching. Such instability has been observed 
and analyzed by Oei, Goldberg and Chang [68]. The remedy recommended in that 


25The Cambrian explosion and various pre-historic mass extinctions notwithstanding [39]. 
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report, termed tournament selection with continuously updated sharing, is very similar 
in principle to the notion of a steady-state process. 

The replacement mechanism chosen here utilizes the three selection sub-routines 
already in place. For the two ranking methods (Goldberg and Fonseca-Fleming), 
the two individuals to be replaced are simply chosen by roulette wheel selection 
based on inverse ranks. For the tournament method, the losers of two independent 
dominance tournaments are chosen for replacement (ties are again resolved by niche 
count). Unlike the parent selection process, replacement does not require the two 
selected individuals to be of the same species. Such stochastic selection of lethals 
should alleviate some of the high selective pressure associated with GENITOR-type 
algorithms. 

The use of steady-state processing combined with clone rejection makes the prob- 
ability of crossover term in the schema theorem (4.1) inapplicable here. In typical 
generational GAs, parents are selected one at a time (by roulette wheel, say) and 
temporarily placed in a mating pool until all selections are complete. It is therefore 
likely that highly fit individuals will be represented several times in the mating pool. 
Those in the mating pool are then paired off randomly, and mated with probability 
pe (the pre-defined crossover probability). If crossover does not occur, the parents 
are simply copied into the next generation after being subjected to mutation. In the 
steady-state process used here, every mating involves crossover and the parents then 
compete with all other members of the population to avoid being replaced by the 
children. All parents are in a sense copied to the following generation, so whether 
this is equivalent to a pe of zero or unity is a matter of interpretation but essentially 


irrelevant. 


4.4 Competitive Species 


The idea of allowing multiple, non-interbreeding species to compete in a genetic algo- 
rithm is an original component of this research. It is intended to allow simultaneous 
optimization of independent and mutually exclusive options. This is accomplished by 
simply placing an equal number of representatives from each species (1.e., each op- 
tion) in the initial population and preventing subsequent inter-species mating. Each 
species will then multiply or diminish relative to the others based on how well it can 
adapt to optimality as defined by the selection method and the evaluator. 

The need for such a mechanism in the current research is due to the incompatibil- 


ity of the four propulsor configurations considered. During the initialization process, 
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the four configurations—rotor-only (one stage), rotor-stator (two stages), stator-rotor 
(two stages), or stator-rotor-stator (three stages)—are all equally likely, and all ac- 
tive propellers must have their circulation distribution specified. Each specification 
requires two inputs—a value at the hub and a value at the tip. Thus the two-stage 
configurations require two more input values than the single-stage configuration, and 
the three-stage configuration requires four more. Chromosomes for individuals of dif- 
ferent configurations must therefore be of different lengths, or must contain unused 
place-holders. 

The chromosome length problem alone does not motivate the species concept. As 
mentioned in the discussion on initialization in Section 4.3.1, neither variable-length 
chromosomes nor dummy values cause insurmountable difficulties, although they can 
be troublesome and inefficient. Rather, the primary motivation for non-interbreeding 
species in this research is the prevention of invalid selection. The problem, as illus- 
trated in Figure 4-3, occurs when a parent passes on dummy (latent) genes which 
then become active in the offspring. The parents in this example are of different 


configurations, and both are assumed to have been selected by some valid, stochastic 


configuration gene latent circulation genes 


CH! 
(valid) 


(invalid) 





Figure 4-3: Latent genes activated by crossover. 


method based on fitness. The first parent (P1) is of configuration 1 (rotor-only), and 
therefore has arbitrary values of hub and tip circulation for its non-existent second 


and third propeller stages. These genes are shaded in the figure, indicating that they 
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are latent, and are marked by an X, indicating that they were ignored by DPLL during 
evaluation and thus had no impact on P1’s fitness value or its subsequent selection. 
The second parent (P2) is of configuration 4 (stator-rotor-stator) and therefore has 
active values of all six circulation genes. If Pl and P2 are mated, and the random 
crossover point falls between the configuration gene and the final circulation gene as 
shown, the result is one child of each configuration. When these children are evalu- 
ated, DPLL will ignore the valid but extraneous circulation genes which the rotor-only 
child (CH1) has inherited from P2. This is, at worst, a bit inefficient; CH1 is a valid 
mixing of the parents’ parameters. What makes the mating invalid is the passing of 
latent, arbitrary genes from P1 to CH2 where they then become active, since CH2 is 
of configuration 4. The problem is worse if the cut location happens to fall within 
the configuration gene itself, as the children will then be of different configurations 
than either parent.”° 

This problem is resolved by the multiple species concept, which is intended to allow 
the four incompatible configurations to be optimized simultaneously in less time than 
would be required for separate optimizations. A first parent is selected from the 
general population without regard to species, using some standard stochastic method 
based on fitness. Selection of its mate, however, is restricted to members of the 
same species. The same selection method is used for the second parent, but only 
individuals of the same species as the first parent are eligible—the potential mating 
pool is segregated. A minor difficulty with this process, namely what to do if only one 
member of a species remains in the population and happens to be selected, is easily 
overcome by implementing an “endangered species” strategy. This simply spares the 
last two representatives of any species from extinction; however, it is probably not 
crucial to overall results. Any species which has been reduced to only two remaining 
representatives is obviously not well-adapted and is unlikely to play a role in further 
convergence of the algorithm. 

Certainly the four propulsor configurations considered here could be optimized 
independently using custom chromosome lengths and a modified algorithm for each. 
This would result in a unique Pareto frontier for each configuration. These frontiers 
could then be overlaid, if desired, to obtain a “meta-frontier” having configuration 
as a location-dependent characteristic rather than a fixed, global parameter (see Fig- 
ure 4-4). This “frontier of frontiers” would present all the information needed—in 
condensed form—for making a final decision, assuming that no external reasons exist 


261t is certainly possible that, by accident, such a child would turn out to be highly fit. On average, however, it 
would be inefficient to spend time evaluating such “illegitimates.” Injecting arbitrary bits into the selection process 
is the job of the mutation operator. 
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Figure 4-4: Notional Pareto meta-frontier, with objectives A and B to be minimized. 


for preferring one configuration over another. 

Production of the meta-frontier while simultaneously solving the latent gene prob- 
lem is the motivation behind the competitive species concept; elimination of inferior 
species from consideration as the population evolves should make this a more efh- 
cient approach overall than performing separate optimizations. The concept may be 
generalizable to include any design or decision where mutually exclusive and GA- 


incompatible options exist. 
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Chapter 5 


Results 


5.1 Comparison of Selection Methods 


Three versions of the Pareto GA described throughout Chapter 4 are run to the 
completion of 3500 cost function evaluations (CFE). These versions differ only in 
the selection method they employ; the selection methods themselves are described in 
Section 4.3.3, page 107. The cost function is DPLL v2.0; the Pareto objectives are 
maximum internal volume and minimum power coefficient for ducted propulsor sub- 
marines. Internal volume corresponds to stern fullness, as all variants are of identical 
length and diameter and use a common bow profile. All results documented here, 
including the three-objective results of Section 5.3, are obtained using a population 


size of 200 and all are started from the same randomly generated population. 


5.1.1 Final Populations 


Figures 5-1, 5-2 and 5-3 show objective space plots of the final populations produced 
by the three versions. The horizontal axis is normalized internal volume (as defined in 
Appendix B) converted to a minimization. Greater volume (increasing stern fullness) 
is toward the origin. The vertical axis is power coefficient, with a smaller value 
indicating that less power is required to maintain the given forward speed. Variants 
identified as infeasible in the plots are those for which flow separation occurs forward 
of the rotor. 

A prominent feature of these plots is the tendency of the solutions to align vertically 
at discrete values of volume. This is a result of parameter discretization and is not 
a characteristic of the GA. Due to the way this particular optimization is set up, 
the volume objective happens to be strongly dependent on one input parameter— 
fullness factor—and weakly dependent on all others (via the gear and shaft volume 


penalties described in Appendix B). Since fullness factor, like all input parameters, 
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Figure 5-1: Initial and final populations from Goldberg ranking method. 
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Figure 5-2: Initial and final populations from Fonseca-Fleming ranking method. 
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Figure 5-3: Initial and final populations from tournament selection method. 


is necessarily discretized, the evaluated volumes are also somewhat discretized. This 
observation prompted the higher resolution of fullness factor relative to the other 
inputs as mentioned in the introduction to Section 4.3.1. 

It is also apparent that DPLL’s power-volume frontier is constrained by the max- 
imum stern fullness allowed in the input. Solutions lying on the constant-volume 
border on the left side of the plots are “weakly” Pareto optimal—they all have nearly 
the same value of one objective. Weak crowding is not a desirable condition in a 
Pareto GA, as it reduces population diversity without supplying any additional infor- 
mation.’ Niching penalties alone, at least as they are used in this research, do little 
to alleviate weak crowding; apparently, the non-dominated status of these solutions 
overcomes their penalties due to niching and they remain likely candidates for selec- 
tion. Tolerances on objective values during dominance checking, intended to make 
some weakly optimal solutions appear dominated, are incorporated in the algorithm 
but are evidently unsuccessful at reducing the effect. 

Generally, a frontier constrained by parameter ranges would indicate that these 
ranges should be increased or adjusted. In retrospect, it is possible that increasing 


' This observation is the author’s own; none of the Pareto GA papers reviewed mention weak crowding. 
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the upper limit on fullness factor in this research would have reduced the weak op- 
timality stacking effect and provided a more complete frontier. On the other hand, 
the maximum fullness allowed here is already quite extreme compared to operational 
designs (see the outer envelope in Figure 4-2) and, since it appears that the minimum 
power coefficient increases rapidly with stern fullness when 1 — Cy < 0.14, sterns of 
greater fullness than the maximum allowed here would probably be very inefficient. 
Comparison of the final population plots from the three methods reveals some 
disagreement regarding the form and location of the non-dominated frontier. Both of 
the ranking methods produce a frontier of very limited span consisting exclusively of 
high-volume variants. The Fonseca-Fleming method appears to have located variants 
of greater efficiency than the Goldberg method, while the Goldberg method retains 
greater volume-diversity in its final population. The general location of the frontier, 
however, is the same for both. Tournament selection maintains the best volume- 
diversity of the three, but this is probably a consequence of the fact that it fails to 
locate the global power minimum at 1 — Cy & 0.16, which would otherwise dominate 
most of the objective space. As will be seen when the frontier’s composition is pre- 
sented in Section 5.2, this is an isolated minimum in terms of propulsor configuration, 
making it difficult for a GA to locate (and, incidentally, very unlikely to be located 
by a hill-climbing method, without some improbably accurate prior assumptions). 
The fact that the true frontier—taken hereafter to resemble that produced by 
the Fonseca-Fleming version—seems limited to very full sterns is hydrodynamically 
interesting and will be discussed in Chapter 6, but leaves much of the (presumably 
dominated) feasibility boundary unresolved. The form of the entire boundary, includ- 
ing the dominated regions, would be useful in supporting the validity of the frontier 
and in learning how minimum attainable power varies across the entire range of stern 
fullness. Fortunately, some indication of the dominated boundary’s location is avail- 
able from the cumulative information produced throughout the GA’s iterations. 
Figures 5-4, 5-5, and 5-6 are the results of overlaying all feasible variants generated 
throughout 3500 reproduction and evaluation cycles and tracing the boundaries of the 
resulting scatter in objective space. This procedure reveals the Pareto frontier in each 
case as well as an approximation of the dominated feasibility boundary. Resolution 
of the region to the right of the minimum at 1 — Cy ®& 0.16 is relatively poor for 
both of the ranking methods because it is dominated and the GA allocates most 
of its computation time to the non-dominated frontier. The tournament method, 
however, in failing to locate the minimum, provides corroborating evidence of the 


lower-volume feasibility boundary as it conducts a much more extensive search in 
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Figure 5-5: Feasibility boundary from Fonseca-Fleming ranking method. 
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Figure 5-6: Feasibility boundary from tournament selection method. 


this region. These incomplete but informative results confirm a slight downward 
trend in power coefficient for very slender sterns which is barely noticeable in the 


feasibility plots of the ranking methods. 


5.1.2 Convergence Analysis 


Defining a measure of performance for a Pareto GA is a difficult matter. Simply count- 
ing the number of non-dominated solutions in the population after a given number of 
CFEs is reasonable, but does not account for the diversity or quality of the solutions. 
Even initial, randomly generated populations have several non-dominated members 
if domination is defined only in terms of contemporaries, which it must be since the 
location of the true frontier is never precisely known. The number of Pareto solutions 
found, therefore, should not be used as the sole criterion for judging performance. 
The quality of a Pareto GA’s frontier actually depends on several criteria, including 
accuracy, resolution, range (i.e., how closely the true limits are approached), and 
density distribution. 

There are few suggestions in the literature for quantifying the performance of 


Pareto methods in general. Schott [80] proposed the variance of the Euclidean dis- 
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tance (in objective space) between adjacent frontier solutions as a measure of a Pareto 
GA’s ability to spread the population evenly. This is certainly a valid measure in gen- 
eral, but is not appropriate for this particular problem. Here, the discretized objective 
space combined with a very limited frontier span produces misleading variance values. 
Schott also proposed a “7-point distance measure,” where points on the objective axes 
are pre-selected and the distances from each point to the nearest feasible member of 
the final population are averaged. As with variance, this is a measure more suited for 
continuous frontiers defined by a statistically significant number of solutions, and is 
not well-suited to this particular problem. 

What can be used here to compare the three selection methods is qualitative anal- 
ysis of the final population and feasibility boundary plots presented above, supported 


by some simple quantitative measures. Figure 5-7 shows characteristics of the pop- 
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Figure 5-7: Comparison of frontier dynamics for the three selection methods. Mean volume, mean 
power coefficient, and number of non-dominated solutions found vs. cost function evaluations. All 
quantities shown involve feasible variants only. 
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ulation, as functions of the number of cost function evaluations, for each method. 
The upper plot tracks the mean volume of all feasible individuals in the popula- 
tion, the middle plot tracks mean power, and the lower plot tracks the number of 
non-dominated solutions as the population evolves. 

Both ranking methods show steady trends in mean objective values; they ap- 
pear to have reached a steady state after approximately 2300 CFE. The number of 
non-dominated solutions, as expected, are erratic and provide little means for dis- 
tinguishing among the methods in this application. Tournament selection shows no 
real trend in mean values and its population is undergoing rapid change at the end of 
3500 CFE. This is probably related to the proportion of infeasibles in the population, 
which for the tournament method begins to increase rapidly around 3000 CFE.” The 
number of non-dominated solutions, as with the ranking methods, is erratic. 

Qualitatively, the ranking methods obviously outperform tournament selection. In 
addition to their steadier trends in population characteristics, they both identify an 
important region of the frontier which is missed by tournament selection. Between the 
two ranking methods, Fonseca-Fleming ranking appears to provide a fuller, smoother 
frontier than Goldberg ranking and, most importantly, identifies several feasible vari- 
ants which would dominate those on the Goldberg frontier if the final populations 
from the two methods were overlaid. 

Although Fonseca-Fleming ranking appears to out-perform the other selection 
methods in all respects, some caution is in order when drawing conclusions from 
these results. Due to the considerable processing time required by this evaluator 
(DPLL) and the time limitations on this research, the total number of CFEs for each 
method is quite small compared to most applications in the literature. It is possible, 
for instance, that tournament selection was close to locating the global power mini- 
mum when it was terminated, and that it would have soon settled there and produced 
a “better” frontier than either ranking method. The results of several undocumented 
preliminary runs, however, indicate that the tournament method does seem to have 


more difficulty in locating this inconspicuous optimal region. 


5.2 Composition of the Feasibility Boundary 


As seen in the previous section, the local position and shape of the final feasibility 
boundary, including the dominated region, varies among the three selection methods. 
2It is possible that the dynamic constraint penalty is somehow incompatible with tournament selection as it is 


implemented here. The rapid increase in proportion of infeasibles was noted for several of the most recent revisions 
of the tournament selection subroutine. 


124 


Despite local differences, however, all methods generate very similar solutions—in 
terms of parameters, that is—at equivalent positions along the boundary. Figure 5-8 
shows a boundary in its (assumed) continuous form, based on the combined results 
from Section 5.1.1. Analysis of interim and final populations indicates that the bound- 


ary consists of three zones defined by propulsor configuration. 
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Figure 5-8: Feasibility boundary composition, from combined results of the three selection methods. 


For gently tapered afterbodies with low volume coefficients (1 — Cy > 0.28), the 
rotor-stator configuration is predominant on the boundary—that is, all of the best 
(lowest power) solutions found in this region throughout the GA’s iterations are 
rotor-stator. Circulation values on the blades of both stages are generally high, but 
do not appear to be limited by the input range. Ducts in this region have a very 
strong tendency toward mid-range chord length (* 0.035) and toward the maximum 
negative load allowed (Gauct = —0.04).2 Mean tip radius is around 0.033L, with little 
variation. Duct axial position with respect to the rotor tip is variable. 

For slightly fuller sterns (the middle zone in the figure), the best propulsor ar- 
rangement in terms of minimal power seems to be stator-rotor. Solutions in this 
region of the boundary tend to have stators with maximum allowed hub circulation 
and near-maximal tip circulation. Rotor hub and tip loads are also high, but slightly 
less, on average, than those of the stator. Duct position relative to the rotor tip tends 
strongly toward the most upstream location allowed (i.e., rotor tip at 80% of duct 
chord); however, duct chord length itself is variable in this zone. The mean duct load 
is negative but of lower magnitude (Gguct = —0.02) than in the first zone and the tip 


3 Negative load corresponds to decelerated flow through the duct and a lift component opposing rotor thrust. 
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radius is slightly decreased. This zone extends to the left past the global power mini- 
mum at 1—Cy = 0.16 and therefore makes up a small portion of the non-dominated 
frontier. 

The third zone, where 1 — Cy < 0.14, is composed entirely of stator-rotor-stator 
arrangements. This zone is entirely non-dominated, although much of its span is 
weakly optimal. Circulation values on the upstream stators are generally of mid- 
range, those of the rotor are somewhat higher (similar in magnitude to the rotor 
values in the middle zone), and those of the downstream stator are below mid-range 
and show very little variance. The duct tends to be positioned such that its mid-chord 
point is slightly downstream of the rotor tip. Mean duct chord length is 0.05£, slightly 
above the median value allowed, with very little variance. Duct load shows essentially 
zero variance with a mean of Gauct = 0 (no net load, but possibly inducing velocities 
near the leading edge to prevent flow separation). Tip radius, as in the other zones, 
averages around 0.03L but shows a slight tendency to increase with stern fullness. 


Most of the non-dominated variants in this zone border on separation. 


5.3 Three-objective Pareto Optimization 


Pareto optimization is not limited to two-objective problems; the definition of Pareto 
optimality is applicable to any number of objectives greater than one. T'wo-objective 
spaces happen to be the easiest by far to display and interpret and are therefore seen 
in the majority of examples and test cases. In general, the Pareto frontier is a surface 
of dimension n — 1, where n is the number of objectives considered. 

This section documents the performance of the Pareto GA in defining a two- 
dimensional frontier, with minimal cavitation index (defined in Appendix A) included 
as the third objective. The changes involved in converting the GA from two objectives 
to three are minor; only the selection subroutines and some minor input/output 
features are affected. The starting point is the same initial population used for all 
the two-objective optimizations presented above, with a population size of 200. 

Figure 5-9 shows the resulting two-dimensional Pareto surface after 2900 CFE, 
when 24 non-dominated solutions exist in the population. These appear as the white 
squares in the figure; the surface is interpolated from these points.* Note that the 
axes of this plot are the same as those for the two-objective results, and that the 
locations of the frontier points correspond closely to the feasibility boundary plots of 
secu 5.1.1. 


4This interpolation was done with Tecplot v7.0, using the “kriging” scheme. 


126 


0.4 


0.35 


0.12 0.16 0.20 0.24 0.28 





E 


y 


Figure 5-9: Non-dominated frontier for three objectives, using Fonseca-Fleming ranking method. 
The symbols are the non-dominated points after 2900 CFE. The cavitation surface is interpolated 
from these points. 


Analysis of the final population shows that all of the non-dominated solutions in the 
region 1 — Cy > 0.18 are rotor-stator configurations. These lower-volume variants 
account for 16 of the 24 solutions in this Pareto set, and include the two lowest 
cavitation indices found. Two stator-rotor configurations are present at 1—Cy ~ 0.15 
and have the lowest power coefficients in the set, but not as low as those of the two- 
objective minima presented above. Presumably, further iterations would locate these 
points, as any solution which is non-dominated in two objectives must remain so 
when a third is added. For larger values of volume, stator-rotor-stator variants again 
predominate; however, at the extreme upper limit of volume are two solutions having 
a rotor-only configuration. In general, lower cavitation indices appear to correlate 
to rotor-stator arrangements and higher indices to stator-rotor-stator arrangements. 
Tip radii throughout the set are below mid-range (œ~ 0.025L) and blade loading is 
relatively light compared to that of the two-objective set. Duct load shows a solid 
trend from slight negative for low-volume variants to mid-range positive (G = 0.05) 


for the fullest sterns with rotor-only propulsors. 
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Chapter 6 


Conclusions 


6.1 Hydrodynamic Issues 


When drawing hydrodynamic conclusions from the results of Chapter 5, it must be 
kept in mind that what is being optimized here is the DPLL model. While much effort 
is made in this research to ensure and verify the revised DPLL’s accuracy, it remains 
a preliminary, exploratory code and is subject to its assumptions and simplifications 


as put forth in Chapter 2. 


6.1.1 Objective Relationships 


The two-objective results presented in Chapter 5 strengthen the case for full stern 
submarines and support the first thesis of this research. Optimization of the DPLL 
model indicates that the power required to propel a submarine at steady speed does 
not increase with stern fullness if the propulsor configuration is matched to the hull 
profile. In fact, it may be that less power is required for fuller sterns, and there may 
be an optimal stern profile, quite full, which minimizes required shaft power. 

The three-objective results indicate that the likelihood and severity of cavitation 
generally increases with stern fullness and propulsive efficiency. Based on interpola- 
tion of non-dominated solutions, it appears that a power trade-off is a more efficient 
means of reducing cavitation than a volume trade-off. 

It may be possible in general to extract frontiers of lower dimension from surfaces 
of higher dimension; for example, the three-objective Pareto surface of Section 5.3 
might allow determination of the power vs. cavitation frontier. This could be done by 
sorting the final non-dominated solutions in terms of two objectives, without regard 
to the third, and extrapolating the results. It is unclear how the accuracy of such 
a process would compare to that of an explicit, two-objective optimization. This 


potentially useful aspect of Pareto optimization is not investigated here. 
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6.1.2 Frontier Composition 


The most interesting features of the frontier and its extension, the dominated feasibil- 
ity boundary, are the changes in propeller arrangement and duct loading with stern 
fullness. Following are possible hydrodynamic explanations for the trends noted. 

Flow near the stern should become more radial as the stern becomes fuller and the 
velocity magnitude should decrease due to the increasing effect of the wake fraction. 
Radial contraction intensifies any tangential velocity present; this is presumably the 
reason that compound propellers, which can take advantage of this effect, make up 
the entire frontier. The decelerated flow (wake fraction) caused by viscous effects 
and negative duct load should increase the rotor’s efficiency, but also makes sepa- 
ration of the hull boundary layer more likely. The trend from negative duct loads 
for low-volume variants to zero load for full sterns evidently results from a tradeoff 
between these two effects. Since separation is not a problem for the low-volume, gen- 
tly tapering sterns, a relatively large negative duct load can be used to reduce the 
rotor inflow velocity. As stern fullness increases, however, negative loads of the same 
magnitude result in separation and are therefore not feasible. This loss of potential 
rotor efficiency is compensated by moving the rotor downstream of the stator. In this 
location, it can benefit from the increase in stator pre-swirl provided by radial con- 
traction. The trend toward zero duct load continues as the stern becomes very full, 
and zero load predominates among Pareto solutions at the upper limit of volume. In 
this zone, any attempt to decelerate the rotor inflow is likely to result in separation. 
This further loss of potential efficiency is compensated by an additional stator. Being 
limited in the magnitude of its circulation values, the rotor is unable to utilize all of 
the tangential velocity imparted by the forward stator and magnified by the radial 
contraction. The downstream, lightly loaded stator probably captures enough of this 
otherwise wasted energy to overcome its viscous penalty. 

It should be noted that the high-volume Pareto solutions in all of these results 
are generally very near separation. In moving from an exploratory optimization such 
as this research provides to more detailed analysis, some margin of safety should be 
implemented to exclude borderline solutions which might separate when such things 


as surface roughness, appendage effects, and unsteady forces are accounted for. 


6.1.3 Power Density on the Frontier 


The two-objective results indicate that for given length and forward speed, a fuller 
stern (and therefore greater volume) does not require additional shaft power and 


may in fact require less, depending on the propulsor configuration. In practice, this 
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advantage would likely be used to shorten the overall length of the submarine while 
maintaining the same volume, rather than simply making the stern “roomier” inside. 
Warren has investigated such shortening and found it to be practical in terms of 
arrangements [88]. 

Figure 6-1 shows the final population from the two-dimensional Fonseca-Fleming 


ranking method overlaid with lines of constant power per unit volume. Clearly, min- 
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Figure 6-1: Initial and final feasible populations from Fonseca-Fleming ranking method, overlaid 
with lines of constant power density. 


imum “power density”—or greatest volume propelled per unit power—of all feasible 
solutions is located near the global Pareto power minimum at 1 — Cy = 0.16. It is 
tempting to conclude that this point represents the optimum, but note that minimiza- 
tion of power density is simply a scalarization and should not be used independently 
of a Pareto method. In this particular case, minimum power density simply provides 
a logical means of choosing from among the solutions on the frontier and eliminates 


the weakly optimal portion of the frontier from consideration. 
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6.2 Optimization and Decision-Making Issues 
6.2.1 General Comments 


The second thesis proposed in Chapter 1—that knowledge of the solution space is 
necessary for informed decision-making—is also validated; the dominance of fuller 
sterns over those which are more tapered is certainly non-intuitive. Such configura- 
tions would probably not be investigated without some indication of their potential, 
as is provided by the frontier. 

It is notable that the ratio of cost function evaluations to design space size used here 
is quite small relative to most published applications. The two ranking methods are 
reasonably stable after 3500 CFE, which is less than one ten-millionth of the design 
space possibilities. This illustrates the power of implicit parallelism in exploring 
objective space and exploiting known solutions. 

The process of defining the dominated feasibility boundary by the cumulative 
results of the Pareto GA is somewhat cumbersome and imprecise, but does provide 
supporting data for the final frontier at no additional computational cost. Definition 
of this boundary may be useful in general Pareto optimizations where the frontier 
itself is quite limited, as it is here. 

The full stern power minimum discovered by the two ranking methods is isolated 
in terms of propulsor configuration. The majority of the Pareto frontier is composed 
of stator-rotor-stator variants; only a small portion, at the lower limit of power coef- 
ficient, is made up of stator-rotor variants. Such a point would be very difficult for a 
gradient-type method to locate, particularly if the four different configurations were 
optimized simultaneously as they are here. 

The three-objective frontier produced by the GA shows how the Pareto concept 
can be extended into as many dimensions as necessary. The results indicate that 
the objective of minimal cavitation is incompatible with both efficiency and stern 
fullness in terms of the DPLL model. In fact, cavitation index is at a maximum in 
the region of the two-objective Pareto frontier. This is likely due to the predominance 
of high rotor blade loads on the two-objective frontier, which tends to result in high 
cavitation indices according to Equation A.1. 

The algorithm would be perfectly capable of evaluating four or more Pareto objec- 
tives, but problems would be encountered in displaying and interpreting results. Prob- 
lems of presentation could possibly be resolved by resorting to “slices” through the 
frontier where one or more objective values are held constant. For several objectives, 


however, interpretation becomes the primary obstacle. Each additional objective 
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provides the members of the population with another way of being non-dominated. 
Eventually the Pareto set becomes so large and diverse that it is of little use to a 
decision-maker. Trends and knees, two of the most useful results of Pareto optimiza- 
tion, become lost in the information overload. This research makes no particular 
attempt to resolve the problems of numerous objectives, other that to propose that 
in some cases the combining of low-level objectives into two, three or possibly four 
meta-objectives may be necessary. Analysis of the results may then indicate better 


ways of combining the objectives to expose trends and knees. 


6.2.2 Selection in Pareto Genetic Algorithms 


Based on the limited results of Chapter 5, the Fonseca-Fleming ranking method 
appears to give the best performance in this application. This conclusion is supported 
both quantitatively by the population’s trend is objective values and qualitatively by 
inspection of the final populations as plotted in objective space. Tournament selection 
is somewhat erratic and is apparently prone to a runaway proportion of infeasibles, 
despite the dynamic penalty threshold implemented here. 

The Pareto GA and its selection subroutines developed during this research went 
through two years of revisions and corrections. The results shown in Chapter 5 are 
from the most recent versions, which were not completed in time to allow numerous 
independent runs, with varying random number seeds, for each selection method. 
Average results over many such runs would be necessary to draw firm conclusions 
regarding the relative performance of the selection methods. However, the final results 
presented are very typical of those seen throughout the development of the GA. In 
other words, despite the limited performance sampling presented here, there is a good 
deal of cumulative, preliminary evidence that the Fonseca-Fleming ranking method 
outperforms the other two, and outperforms tournament selection by a considerable 


degree. 


6.2.3 Competitive Species 


The competitive species concept originated here is very promising. Its application 
reveals that the optimal number, type and relative positioning of propeller stages are 
dependent on stern fullness; the meta-frontier of the DPLL model has in fact been 
identified. Undoubtedly, compilation of the meta-frontier via separate optimizations 
of all configurations would be less time-efficient, although this assertion is not proven 
by direct comparison in this research. Such a comparison would require separate, 


tailored versions of the GA for operating on the different configurations’ chromosome 
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lengths. This would be an interesting experiment, however, particularly if the results 
were overlaid (e.g., Figure 4-4) for comparison to the meta-frontier obtained here. 
Analysis of the final populations from all three selection methods reveals that each 
propulsor configuration is present in approximate proportion to the percentage of the 
boundary it comprises. This indicates that the configuration zones are valid, insofar 
as the DPLL model is concerned, and that the solutions there are not sub-optimal 


accidents residing on a false boundary due to lack of sufficient exploration. 


6.2.4 Miscellaneous 


Penalties due to niching and constraint violation seem to be correctly applied, judg- 
ing by the composition of the final two-objective populations of Chapter 5. The 
number of infeasible solutions in the population remains fairly constant for the two 
ranking methods, supporting the validity of the dynamic penalty scale described in 
Section 4.3.3, page 106. The tournament method typically experiences a sudden, 
uncontrolled rise in the proportion of infeasibles at some point during the iterations. 
This is probably due to some incompatibility between this type of selection and the 
dynamic constraint penalties applied, but the effect has not been investigated in de- 
tail. Crowding is most evident at the weak Pareto boundary, but is probably due in 
part to the somewhat discrete objective space of this problem. It is not clear that this 
phenomenon can be prevented in the general case without some application-specific 


tailoring. 


6.3 Future Work 


As with most time-limited projects, there are areas here which warrant further effort 
or investigation. These fall into two general categories: improvements and modifica- 
tions to DPLL, and further research into the use of Pareto GAs for optimization and 


decision-making. 


6.3.1 Evolution of DPLL 


Following are possible modifications to DPLL which in the author’s estimation could 


result in greater accuracy, speed, and robustness. 


1. A front-end routine is needed to convert body and duct geometries to B-spline 
vertex files for DPLL input. Modeling an existing body and/or duct in the cur- 


rent version requires tedious trial and error to get the vertices correct. ‘The user 
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should be given the option of providing vertices or offsets; conversion, if neces- 
sary, should take place internally. This will probably require some interaction, 
since the best possible fit using a given number of vertices varies according to 


the complexity of the prescribed curve. 


. The body model should be converted from submerged source rings to a paneled 
surface. In the current version, source ring strengths can be erratic depending 
on the number and spacing used. Converting to a panel method may result in 
better pressure drag calculations, which have been somewhat suspect throughout 


this research. 


. The duct and blade representations should be given thickness distributions, by 
adding sources to the system matrix which solves for singularity strengths. The 


additional computation load should not be significant. 


. The duct boundary layer should be modeled, in a manner similar to that of the 


body, to predict leading edge cavitation or separation. 


. Tip gap and blockage effects should be included. In keeping with the intended 


use of the code, a parametric estimate is probably sufficient. 


. Linear circulation distributions on lifting lines are restrictive. Goldstein calcu- 
lations should be modified to allow specification of individual circulation values 


on all segments of the lifting lines. 


. The body B-spline should update at each iteration so that the radius of the 
sting equals the calculated radius of the hub vortex. This will prevent wake 
streamlines from entering the hub vortex region, a condition which currently 


must be monitored and corrected manually. 


. Submarine propeller blades are usually quite skewed; DPLL’s lifting line model 


should be modified to allow more precise specification of skew and rake. 


. More validation against published data is needed. 


6.3.2 Optimization 


Population size was originally intended to be a variable in this research, but selection 


method was considered more important and there were insufficient resources of time 


and hardware to investigate both issues. It would be interesting to know the GA’s 


convergence rate as a function of population size, and to determine the impact of the 


multiple species on optimal population size. 
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Two parameter ranges—duct load and stern fullness—apparently constrained the 
frontiers in this research. Although the allowed limits on both of these seem to 
already be quite extreme, a non-constrained frontier would be interesting. Properly 
adjusted parameter ranges might eliminate the weak crowding observed in the two- 
objective frontiers. If not, some mechanism should be devised for reducing this effect. 
One possibility is to apply negative dominance tolerances to solutions below the first 
dominance layer, thus increasing their apparent dominance count. 

The results given by this exploratory model should be investigated with more 
sophisticated codes. In particular, the full stern, stator-rotor variants with slight 
negative duct load, which represent the apparent global power minimum, should be 
analyzed in greater detail. 

The Pareto GA used here is not the only possible method for locating non- 
dominated solutions, although it seems to be the only way to locate several in a single 
run. This method should eventually be compared in terms of efficiency and efficacy 
to some baseline, such as stochastic hillclimbing, and thereafter to other alternatives, 


such as systematically varied objective weighting. 
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Appendix A 


Cavitation Index 


Farrell and Billet [25] cite three types of cavitation which occur in the tip region of 
an axial-flow pump: gap, blade-end and leakage vortex cavitation. Gap cavitation 
occurs when flow separation occurs in the gap between the blade tip and the housing or 
casing.’ The likelihood of gap cavitation may be reduced by rounding the intersection 
of the blade surfaces with the tip [29]. Blade-end cavitation is due to boundary layer 
vorticity stretching and occurs when tip clearance is small. DPLL does not model 
boundary layer profiles on the blades or the duct; therefore, this type of cavitation 
cannot be taken into account here. Leakage vortex cavitation is due to low pressure 
within vortices formed by the interaction of through-flow and gap flow. 

According to Farrell and Billet, this is the “most problematic and least understood” 
of the three types. Based on experimental results? and comparison with published 
data sets, they suggest a correlation between relevant design parameters and the 
minimum leakage vortex pressure. This correlation is adapted here to predict relative 
likelihood of cavitation among variants in the optimization process: 

Pike Wi (1 — e7***) + Wook, i 

8m? A? (1 — en U 
The correlation parameters are defined in Table A.1, along with their adaptations in 


the DPLL model. The usual definition of cavitation index involves the total far-field 


on Re! (A.1) 


o 


pressure pœ, the vapor pressure at the operating temperature p,, and the far-field 
velocity (assumed here, without loss of generality, to be ship speed V;): 
aie (A.2) 
PVs 
This is an environment-dependent positive quantity which decreases with increasing 
reference velocity or decreasing po. A lower value of ø indicates greater likelihood of 
Ee a enone cle to ducted propellers. In adapting the Farrell and Billet correlation 


to this research, the housing is considered to be the inner surface of a duct. 
?These experiments were performed at the HIREP facility discussed in Section 3.1.2. 
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Parameter Farrell and Billet DPLL 
K correction factor for levels of air con- wo 
tent in water 
A empirical constant (0.36) 0.36 
W, relative velocity at blade tip relative velocity at control point of 
outer lifting segment 
Wao edge velocity of end-wall boundary Wi 
layer at blade tip 
(Che tip lift coefficient at zero clearance axial force produced by outer lifting 
condition segment normalized by ¿pW A 
€ ratio of tip gap to blade chord length 1% of tip radius divided by blade 
at tip chord length, which is constant for 
all variants at 2% of body length 
ES empirical constant (0.18) 0.18 
A ratio of tip gap to maximum blade 1.0 
thickness at tip 
U tip speed propeller angular velocity x radius 
of outer lifting segment control point 
Re, blade tip chord relative Reynolds Weo = Wj, blade chord length con- 


stant at 2% of body length, viscosity 
derived from input vehicle Reynolds 
number “£ 


number, o 





Table A.1: Adaptation of cavitation parameters to DPLL model 


cavitation. The local pressure coefficient is defined somewhat similarly as the differ- 
ence between local total pressure and reference pressure, normalized by a convenient 
dynamic pressure: 
P — Poo 
e == ee 
PVs 


(A.3) 


Obviously, local cavitation is dependent on local pressures. If the local pressure at 
some point in the flow field is less than vapor pressure, cavitation will occur. On 
the other hand, if the minimum pressure of the entire flow field is greater than 
vapor pressure, cavitation is unlikely. If this is the case, and if the global cavitation 
index is then somehow steadily lowered, cavitation will first appear in the flow when 
o= C, 
for which cavitation will occur. When comparing propulsor configurations operating 


(that 


is the more desirable in terms of cavitation, as it requires a 


In other words, the quantity —C,_,, is the maximum cavitation index 


min’ min 


in the same global conditions, therefore, the one having a higher value of Cp, 
is, a lower value of —Cpmin 
lower global cavitation index (i.e., lower static pressure or greater speed) to perform 
as poorly. This is why the ø of Equation (A.1), as calculated in the DPLL model, is 
to be minimized if minimal cavitation is desired. 


In this research, it is assumed that a direct correlation exists between the calcu- 
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lated cavitation index and propeller noise level. There are several reasons, however, 
why this must be treated only as a relative, first order approximation. First, some 
parameters in Equation (A.1) are not calculated by the DPLL model and their values 
must either be ignored or assumed in the adaptation, as shown in Table A.1. Also, as 
mentioned above, the leakage vortex cavitation estimated by this correlation is only 
one of at least three types known to occur in the tip region. No account is taken of 
these other cavitation mechanisms. Finally, cavitation is certainly possible in regions 


other than the tip and the correlation makes no prediction in this regard. 
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Appendix B 


Reduction Gear and Shafting 
Volumes 


Reduction gear volume is a function of the diameter and face width of the gears. The 


maximum tangential load W; per unit length on the gear teeth can be written as [42]: 


(FE)ma =a ide (B.1) 
where «x is a safety factor, J is an experimentally determined material constant in 
units of force per unit length of face per unit length of diameter, d, is the gear 
diameter, and F, is the effective face width of the gear. The actual working load is a 
function of power and shaft speed: 

W: Z 
F, wd,F. 


where P is power delivered by the prime mover and w is the angular velocity of 





(B.2) 


the gear. According to [42], the most “economical” reduction gear uses the smallest 
pinion diameter possible in relation to its working face. However, the face width-to- 
diameter ratio must be small enough to avoid excessive deflections; a typical value is 
2.25. In any regard, the face width is of the same order as the diameter, so by setting 


the two expressions equal one may write gear diameter in terms of torque (Q): 
OS Qs (B.3) 
and the projected area of the gear is then 
As = kQ? (B.4) 


where k, is constant for a given gear material, width-to-diameter ratio and safety 


factor. 
It is necessary to approximate k, in this research so that valid volume penalties 


may be applied to variants with high torque coefficients. To avoid the inclusion of 
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classified information, a published study of general submarine internal arrangements 
is used as a baseline [16]. The mechanical drive variant described there produces 
25,000 hp at 100 rpm, resulting in a torque of 1.79 x 10% N- m. The cross-sectional 
area of the reduction gear for this variant, assumed to be representative of those 
currently in use, is approximately 50 m?, giving a k, of 0.0034. This constant, along 
with the assumption of a length-to-diameter ratio, allows approximation of reduction 
gear volume given torque. 

The volume required for the shaft is calculated similarly, although shaft length is 
assumed constant for all variants and is thus included in the constant of proportion- 


ality. This reduces the torque contribution by an order of magnitude: 
A, = k,Q? (B.5) 


Using arrangements drawings from [16] to estimate the shaft cross-sectional area 
and the powering parameters from above, k, is found to be 0.075. 

The volume coefficient Cy used in this research is defined as the volume enclosed 
by the hull profile minus the estimated gear and shaft volumes, normalized to a right 


circular cylinder of radius Ag and length 1.0. 
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